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I 16.  Abstract  (Concl) 
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] large  scale  simulations  of  actual  crash  tests. 
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1 . INTRODUCTION 


The  increased  concern  in  recent  years  over  vehicle  safety  has 
focused,  in  part,  on  the  ability  of  a vehicle  to  sustain  a crash 
event;  to  absorb,  redirect  and  otherwise  manage  the  severe  forces 
and  energy  associated  with  the  event;  and  thereby,  to  lessen  the 
severity  of  the  environment  to  which  passengers  are  exposed.  In 
addition  to  basic  vehicle  crashworthiness,  serious  attention  has 
also  been  devoted  to  such  factors  as  the  propensity  of  a vehicle 
to  damage  in  minor  crash  events  and  the  threat  posed  to  pedestrians 
or  other  vehicles  by  particular  vehicle  designs. 

Since  quantitative  knowledge  of  the  deformation  characteris- 
tics of  vehicle  structures  is  an  intrinsic  requirement  in  all  con- 
siderations of  this  type,  it  is  not  surprising  that  the  developing 
concern  for  vehicle  safety  has  given  rise  to  intense  and  continuing 
efforts  to  provide  mathematical  descriptions  of  structural  deforma- 
tions in  the  crash  environment.  This  report  presents  the  results 
of  a research  project  undertaken  by  I IT  Research  Institute  (IITRI) 
for  the  Department  of  Transportation  to  develop  a finite  element 
computer  program  for  use  in  the  dynamic  analysis  of  vehicle  struc- 
ture, including  sheet  metal,  in  a crash  environment. 

1.1  Previous  Analyses 

The  literature  relating  to  vehicle  crash  simulations  is  exten- 
sive and  has  been  summarized  by  several  previous  investigators 
(Refs.  1 through  5).  Present  analytical  procedures  fall  largely 
into  two  broad  categories  which  may  be  termed  lumped  parameter  or 
equivalent  systems  and  formal  approximation  techniques  (primarily 
finite  element).  In  the  former  category,  the  vehicle  is  treated 
as  an  assemblage  of  lumped  masses  and  resistance  elements  which  are 
selected  to  represent  specific  subassemblies  and  components.  The 
properties  of  the  various  elements  of  such  models  are  determined  by 
a static  and  dynamic  crush  testing,  direct  measurement,  and  the  ju- 
dicious use  of  simple  analysis  procedures.  The  use  of  such  tech- 
niques has  become  widespread  in  studies  of  vehicle  crashworthiness 
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as  represented  by  Tani  and  Emory  (Ref.  6),  Kamal  (Ref.  7),  and 
the  models  developed  by  Battelle  (Ref.  8)  and  Dynamic  Sciences 
(Ref.  9).  Such  techniques  are  valuable  in  that  they  provide  an 
efficient  and  inexpensive  method  of  simulation,  are  fairly  accu- 
rate for  well  established  vehicle  configurations  and,  in  effect, 
serve  as  a depository  for  previous  test  and  design  results.  The 
limitations  of  this  type  of  procedure  are  principally  that  they 
are  heavily  reliant  on  previous  testing;  largely  limited  to  well- 
established  vehicle  configurations  and  only  weakly  tied  to  basic 
principles  of  structural  behavior. 

The  finite  element  approach  represents  an  attempt  to  make 
up  for  the  limitations  inherent  in  the  lumped  parameter  analyses, 
by  employing  more  formal  approximation  techniques  in  the  discre- 
tization of  the  structure  and  a greater  reliance  on  more  fundamental 
structural  representations  and  properties.  The  limitations  of  this 
approach  are  found  in  the  inherent  tendency  toward  more  complicated 
and  expensive  computations  and  the  difficulties  found  in  modeling 
the  extensively  complex  phenomena  associated  with  the  crash  envi- 
ronment. Such  phenomena  include  large  deflections  and  rotations 
in  the  deformed  structure,  regions  of  intense  curvature  (wrinkling), 
material  strain  rate  effects  and  interference  and  contact  among 
structural  components  during  the  response.  These  procedures,  in- 
deed, are  not  totally  free  of  a reliance  on  testing  and  analytical 
judgment  and,  in  fact,  can  be  viewed  as  a somewhat  more  formal 
lumped  parameter  system  and  used  in  connection  with  the  more  empiri- 
cal procedure  in  hybrid,  combination  models. 

A variety  of  finite  element  analyses  directed  toward  the  dy- 
namic analysis  of  vehicle  structures  have  been  reported  previously. 
Selna  and  Salinas  (Ref.  10)  report  a three-dimensional,  linear  elas- 
tic dynamic  analysis  of  a two-door  sedan  in  a 60  mph  frontal  bar- 
rier impact.  Their  results  indicate  that  significant  material 
yielding  and  collapse  would  occur  and  clearly  demonstrate  the  body 
and  engine  pitch  motions  which  ensue  in  such  a crash.  Kirioka  (Refs. 
11  through  13)  in  a series  of  papers  have  reported  an  elaborate  com- 
putational system  for  the  analysis  of  vehicle  structures.  The 
analysis  includes  nonlinear  material  and  geometric  effects,  and 
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static,  dynamic,  and  modal  procedures  and  employs  graphic  sys- 
tems and  mesh  generation  as  analysis  aids.  Thompson  (Ref.  14) 
has  investigated  the  side-on  collision  of  passenger  vehicles 
using  finite  element  beam  elements  to  model  the  vehicle  structure. 
The  elements  include  nonlinear  material  properties  and  a geometric 
stiffness  matrix  for  beam  columns  and  shows  reasonable  correlation 
with  test  data  despite  the  lack  of  plate  and  membrane  elements  in 
the  analysis.  Shieh  (Ref.  15)  has  developed  a framework  analysis 
for  the  precollapse  dynamic  response  of  plane,  ideal,  elastic- 
plastic  frame  structures  which  he  applies  to  the  analysis  of  frame 
work  devices  used  to  attenuate  crash  effects.  Finally,  Young 
(Ref.  16)  and  Melosh  (Ref.  17)  have  reported  on  an  analysis  employ 
ing  an  energy  minimization  technique  for  large  displacements  of 
elastic-plastic  beam  and  truss  elements.  Results  reported  in 
Ref.  17  suggest  that  although  the  equivalent  framework  used  is 
substantially  stiffer  than  the  actual  structure  the  computed  re- 
sponses are  similar  in  character  to  those  obtained  in  test. 

1.2  Analysis  Requirements 

The  response  of  vehicle  structures  under  crash  loadings  is  a 
complex  process  primarily  involving: 

« transient,  dynamic  behavior 

0 complicated  framework  and  shell  assemblages 

0 large  deflections  and  rotations 

0 extensive  plastic  deformation. 

Previous  attempts  at  a formal  analysis  of  this  process  have  been 
only  partly  successful  due  to  a variety  of  limitations  which,  in 
particular  instances,  have  included  inadequacies  in  element  formu- 
lations, material  representations  or  solution  procedures.  The 
work  presented  in  this  report  represents  an  attempt  to  develop  a 
finite  element  program  which  is  specifically  tailored  to  the  class 
of  problems  inherent  in  vehicle  crash  response,  and  which  employs 
or  extends  current  avenues  in  finite  element  analysis  which  seem 
best  suited  to  such  problems.  The  field  of  nonlinear  finite  ele- 
ment analysis  is  currently  an  extremely  active  area  of  research 
with  an  extensive,  related  literature  and  a variety  of  methods 
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and  approaches.  Consequently,  a formal  review  of  the  field  as 
background  for  the  present  analysis  approach  is  not  attempted. 
Instead,  major  features  of  the  present  work  are  briefly  described, 
and  some  rationale  is  offered  for  their  use  in  the  context  of  ve- 
hicle analysis. 

1.2.1  Solution  Procedure  - The  analysis  procedure  is  based 
on  explicit,  timewise  numerical  integration  of  the  equations  of 
motion  for  the  node  points  (three  translations  and  three  rotations 
per  node).  The  choice  of  an  explicit,  rather  than  implicit,  proce- 
dure is  exploited  throughout  the  analysis  by  using  "lumped"  nodal 
masses  and  by  calculating  the  internal  forces  at  the  nodes  from 
direct  integration  of  the  element  stress  fields  without  reference 
to  element  or  assembled  stiffness  matrices  for  the  structure.  The 
relative  merits  of  implicit  and  explicit  procedures  and  the  im- 
portance of  the  mass  formulation  have  been  studied  by  Krieg  and 
Key  (Ref.  18),  and  the  use  of  a direct  calculation  of  internal 
forces  has  been  described  by  Belytschko  and  Chiapetta.  (Ref.  19)  and 
by  Oden  (Ref.  20).  The  choice  of  an  explicit  integration  procedure 
and  a direct  calculation  of  nodal  forces  results  in  a program  with 
minimum  (but  still  substantial)  computer  storage  requirements. 

This,  in  turn,  is  equivalent  to  a capability  of  processing  reason- 
ably detailed  and  extensive  structural  models  with  relative  ease. 

An  early  capability  for  handling  relatively  large  problems  is  con- 
sidered* essential  in  developing  realistic  simulations  of  actual 
crash  events  because  of  the  complex  geometry  and  construction  found 
in  real  vehicles.  Furthermore,  this  procedure  lends  itself  to  a 
reasonably  well  organized  modular  program  which  admits  subsequent 
change  and  development,  including  the  further  incorporation  of  an 
implicit  method. 

1.2.2  Element  Formulat ion  — The  treatment  of  large  displace- 
ments and  rotations  employs  a decomposition  of  the  element  displace- 
ment field  into  a rigid  body  rotation  and  translation  associated 
with  a local  coordinate  system  attached  to  and  moving  with  the  ele- 
ment, and  a remaining  displacement  field  which  describes  the  de- 
formation of  the  element  relative  to  the  current  position  of  the 
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element  axes.  This  transformation,  in  effect,  removes  the  average 
rigid  body  rotation  of  the  element  and  allows  the  use  of  small  or 
moderate  deflection  element  formulations  in  the  calculation  of 
element  forces.  In  this  manner,  extremely  large  rotations  and 
deflections  can  be  accommodated  by  the  analysis  with  accuracy  de- 
pending primarily  on  the  size  of  the  elements  relative  to  the  curva- 
ture of  the  structure.  Although  a formal  convergence  theorem  is 
lacking,  the  decomposition  does  hold  for  infinitesimal  regions  and 
numerical  studies  show  excellent  agreement  with  classical  solutions. 
Previous  applications  of  this  technique  may  be  found  in  work  by 
Argyris,  Kelsey  and  Kamel  (Ref.  21);  Wempner  (Ref.  22);  Murray  and 
Wilson  (Ref.  23);  Belytschko  and  Hsieh  (Ref.  24);  and  other  refer- 
ences which  are  cited  in  Section  2.  The  computer  program  at  pres- 
ent includes  low  order  triangular  plate  elements  and  three- 
dimensional  beam  elements.  A triangular  membrane-hinge  line  ele- 
ment is  also  available  but  is  not  presently  compatible  with  the 
beam  formulation. 

1.2.3  Material  Properties  —The  computer  program  currently  uses 
simple  elastic-plastic  stress  strain  laws;  a uniaxial  relation  for 
beam  elements  and  a biaxial  strain  hardening  Mises  model  for  plates. 
Element  forces  and  bending  moments  for  given  strain  fields  are  cal- 
culated by  piecewise  linear  numerical  integration  of  the  stresses  at 
selected  points  in  the  cross  section.  The  program  is  designed  to 
accommodate  other  constitutive  relations  readily  or  to  accept  ex- 
plicit relations  among  thrust,  moment,  extension  and  curvature. 

In  fact,  linear  relationships  are  provided  at  this  level  which  re- 
sult in  more  efficient  computations  for  this  class  of  problems. 

1.2.4  Program  Validation  — An  extensive  series  of  test  problems 
were  investigated  throughout  the  program  development  to  ensure 
that  the  basic  formulation  as  well  as  the  numerical  procedures  and 
programming  were  performing  properly.  These  range  from  fairly  simple 
cases  such  as  beam  bending,  wave  propagation  and  arch  buckling  prob- 
lems which  were  designed  to  check  the  basic  formulations,  through 
more  complex  situations  for  which  previously  published  numerical 
solutions  could  be  found.  Finally,  a large  demonstration  was 
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attempted  consisting  of  the  simulation  of  actual  vehicle-barrier 
tests.  In  a subsequent  section,  the  test  cases  are  summarized 
and  the  large-scale  simulation  is  described  in  detail. 

1.3  Report  Organization 

Section  2 of  this  report  contains  a description  of  the  an- 
alysis formulation  in  fairly  complete  detail.  In  Section  3 the 
test  problems  are  summarized  and  a detailed  description  of  the 
large-scale  vehicle-barrier  test  simulation  is  provided.  A sum- 
mary of  the  results  of  the  investigation  and  recommendations  for 
further  application  and  development  of  the  computer  program  are 
given  in  Section  4.  Appendix  A gives  the  complete  plate  element 
curvature  functions  and  detailed  descriptions  of  the  program  input, 
program  operations  and  program  listings  are  provided  in  Appendixes 
B,  C,  and  D. , respectively. 
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2.  ANALYSIS  FORMULATION 


This  description  of  the  analysis  formulation  is  intended  to 
roughly  parallel  the  sequence  of  computation  in  the  computer  pro- 
gram. A typical  program  cycle  is:  integrate  equations  of  motion; 

transform  displacement  to  element  coordinate  systems;  calculate 
element  forces;  and  transform  forces  to  global  and  nodal  coordinate 
systems.  Solution  procedures  employed  in  the  computer  program  are 
also  described. 


2.1  Equations  of  Motion 

The  equations  of  motion  at  a typical  node  point  for  trans- 
lations and  rotations  respectively,  are  written  as  follows: 


Translations 


[M][ug]  = [f|]  - [f*] 


where 

[M] 

[uG] 


the  diagonal,  lumped  mass  matrix  for  the  node 

[ii  ,-,,u  ~,ii  ^,1  is  the  acceleration  vector  for 
L xG  yG  zGJ 

the  node  referred  to  a global  system  [ , Z^,  ] 

common  to  all  nodes 


E 

[ft^]  - the  vector  of  external,  applied  loads  at  the 

VJ 

node  referred  to  the  global  system 
[f^]  - the  vector  of  forces  at  the  node  in  the  global 

Cj 

system  which  results  from  deformations  in  all 
elements  associated  with  the  node. 


(1) 


Unless  otherwise  specified,  all  equations  and  variables  refer  to 
a typical  node  or  element,  depending  on  context,  rather  than  to 
the  ensemble  for  all  nodes  or  elements.  Since  there  is  no  global 
stiffness  matrix  the  formation  of  the  ensembles  is  relatively 
simple. 
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Rotation 
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(2b) 

(2c) 


where  all  components  are  referred  to  a coordinate  system  [x,y,z] 

which  coincides  with  the  principal  axes  of  the  lumped  mass  moments 

of  inertia  at  each  node  and  which  rotates  with  the  node  (i.e.,  a 

set  of  rigid  body  axes  for  each  node)  and  where  [1,1  ,1  ] are  the 

x y z 

principal  mass  moments,  [a  , a , a ] are  instantaneous  angular  ac- 

x y z 

celerations,  [ co  , 03  , 03  ] are  the  angular  velocities,  and  where  [m  ] 
_q  x y z 

and  [m  ] are  the  external  and  internal  moment  vectors  at  the  node 


point 


The  orientation  of  the  nodal  axes  [x,y,z]  with  respect  to  the 
global  axes  at  any  time  during  the  motion  is  established  by  the 
components  of  three  unit  vectors  b^,F^,b^  which  remain  fixed  along 
the  nodal  axes  [x,y,z],  respectively,  as  the  node  rotates.  If  the 
global  components  of  these  three  unit  vectors  are 


(b 

(b 

(b 


— ,b  — ,b  — ) 
xx  yx  zx 

— ,b  — ,b  — ) 
xy  yy  zy 

— ,b  — ,b  -) 
xz  yz  zz' 


(3a) 

(3b) 

(3c) 


then  the  components  of  any  vector,  v,  transform  from  the  nodal  to 
global  system  by  the  transformation 


[vG]  = [B]  [v] 


(4) 


where  the  columns  of  [B]  are  the  components  of  the  three  unit  vectors 


bXK 

b - 
xy 

b — 
xz 

b — 

b - 

b - 

yx 

yy 

yz 

b - 
zx 

b - 
zy 

b - 
zz 

(5) 


Initially,  the  elements  of  the  [B]  matrix  are  obtained  from  the  orienta- 
tion of  the  principal  axes  of  the  mas s moment s of  inertia  for  the  node 
in  the  undeformed  structure  and  are  subsequently  calculated  from  the 
integration  of  the  rotational  equations  of  motion';  equations  (2)  . 
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2.2  Time  Integration 


The  numerical  technique  employed  to  integrate  the  equations 
of  motion  consists  of  the  Newmark-beta  method  (Ref.  25).  This 
method  relates  displacement,  velocity,  and  acceleration  at  the 
beginning  (u0>v0>ao)  and  end  (up  v^.a-^)  of  a time  interval,  h, 
by  the  relations 


where  $ is  an  assumed  parameter  related  to  an  assumed  variation  of 
acceleration  with  time  over  the  interval.  Solutions  to  the  equa- 
tions of  motion,  equations  (1)  and  (2) , may  be  obtained  for  a 
proper  choice  of  3 by  using  equations  (6),  (1),  and  (2)  in  an 

iterative  manner,  predicting  u-^  and  v-^  by  means  of  equations  (6) 
and  correcting  a-^  by  means  of  equations  (1)  and  (2)  until  the  proc- 
ess converges.  The  full  solution  procedure,  with  iteration,  is 
available  in  the  computer  program  but  early  trials  indicated  that 
the  most  efficient  process  was  obtained  with  3=0  and  no  iterations, 
and  all  subsequent  computations  have  been  carried  out  with  these 
parameters . 

The  solution  procedure  is  applied  to  the  equations  of  motion 
for  translational  degrees  of  freedom,  equation  (1),  exactly  in  the 
manner  outlined  above.  The  treatment  of  the  rotational  degrees  of 
freedom,  equations  (2),  differs  slightly  in  that  the  rotation  term 
(analogous  to  uq  or  u-^)  has  no  physical  significance  as  calculated 
and  is  not  subsequently  used  in  the  calculations.  Instead  of  an 
explicit  calculation  of  nodal  displacements  or  Euler  angles,  an 
auxiliary  numerical  integration  is  employed  to  calculate  the  current 
value  of  the  [B]  transformation  based  on  its  previous  value  and  the 
current  values  of  angular  velocity  and  acceleration.  Applying  the 
integration  algorithm  (with  3 = 0)  to  a typical  unit  vector,  b^  , 
results  in 


(6a) 


(6b) 


b 


3 t+h 


= b 


3 t 


(7) 
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The  derivatives  of  the  unit  vector  are  related  to  the  angular 
velocity  and  acceleration  by  the  relations 

(8a) 


b ^ 00  x b ^ 

b^  = oo  x b^+oo  x b^ 

= oo  x + oo  x oo  x 5"^ 


(8b) 


Expanding  into  scaler  components  in  the  previous  [x,y, z] system, 
these  relations  become 

, 2 


b3x 


b,,- 

3y 


= + oo  h + — 7“  (oo  u>  + a ) 
t+h  y 2 x z , y 

h2-- 

. , = - oo  h + -K-  (oo  oo  + a ) 

t+h  x 2 v y z x7 


(9a) 

(9b) 


with  the  third  component,  b^— , obtained  from  the  condition  that 
b^  remains  a unit  vector.  A similar  process  is  applied  to  another 
of  the  unit  vectors,  say  b>2,  and  the  third  vector  is  obtained  from 
the  cross  product  relation,  b^  = b2xb^.  The  components  of  the 
new  unit  vectors  referred  to  the  global  coordinate  system  then 
constitute  the  current  transformation  matrix  [B] . 


2.3  Element  Coordinate  Systems  and  Deformations 

A coordinate  system  [x,y,z]  is  embedded  in  each  element  and 
serves  to  define  the  rigid  body  rotation  of  the  element  and  as  a 
reference  for  element  distortions  and  forces.  The  components  of  a 
vector,  V,  transform  from  the  element  coordinate  system  to  the 
global  coordinate  system  by  the  time  varying  transformation 

[vG]  = [E] [v]  (10) 

where  the  elements  of  E are 


e ~ 

e ~ 

e - 

XX 

xy 

xz 

e ~ 

e ~ 

e ~ 

yx 

yy 

yz 
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e ~ 

e ~ 

zx 

zy 

zz 
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which  correspond  columnwise  with  the  global  components  of  unit 


vectors 

el’e2,e3  oriented  with  the  element  axes 

[£,y,z]  respec- 

tively , 

that  is 



el: 

[ e ~ , e ~ , e ~ J 
xx  yx  zx 

(12a) 

e2: 

[ e ~ , e ~ , e ~ ] 
xy  yy  zy 

(12b) 

e3: 

[ e ~ , e ~ , e ~ ] 
1 xz  yz  zz 

(12c) 

2.3.1  Plate  Elements  —The  element  coordinate  system  for  the 
three-dimensional  plate  element  (Figure  1)  is  constructed  from 
the  current  positions  (i.e.,  undeformed  position  plus  displacement) 
of  the  three  points  (I,J,K)  which  define  the  element.  In  terms 
of  unit  vectors  pointing  from  I to  J (S^j)  and  from  I to  K (SIK) 
the  global  components  of  the  three  unit  vectors,  ei’e2,e3  are 
calculated  as 


el  ^SIJ  + SIK^2 

^ 

e3  = SIJ  x SIK 
e^  = ^ a ]_ 


(13a) 

(13b) 

(13c) 


These  operations  apply  equally  well  to  both  the  deformed  and  un- 
deformed positions  of  the  element  and,  in  particular,  define  e^, 
the  normal  to  the  undeformed  plate  surface,  which  is  subsequently 
employed  to  establish  element  rotations  in  the  deformed  position. 

Element  displacements  and  rotations  are  defined  with  respect 
to  the  element  coordinate  system  and  are  determined  from  the  global 
components  of  displacements  and  the  transformations  [B]  associated 
with  the  three  node  points  I,J,K.  Having  determined  the  [E]  trans- 
formation as  described  above  the  global  components  of  displacement, 
Uq,  at  each  of  the  nodes  as  determined  from  the  integration  of 
equation  (1)  can  be  transformed  to  the  element  system  by  the  oper- 
ation 

[u]  = [E]T[ug]  (14) 

where  [E]  denotes  the  transpose  of  [E] . Since  the  three  node 
points  uniquely  determine  the  transformation  [E]  this  calculation 
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Three-Dimensional  Plate  Element 
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would  eliminate  the  component  of  displacement  normal  to  the  plane 

of  the  plate,  u =0.  For  a small  deflection  plate  formulation 

then  the  nodal  displacements  u , u lie  in  the  plane  IJK  and  de- 

x y 

termine  the  membrane  midsurface  strains  while  the  plate  curva- 
tures are  determined  by  the  rotations  alone.  The  computer  pro- 
gram actually  employs  a somewhat  different,  but  equivalent  and 
more  efficient,  approach  to  the  calculation  of  membrane  strains. 
The  lengths  of  the  element  sides  in  the  deformed  position  are 
calculated  directly  from  the  global  displacements  and  these  new 
lengths  determine  the  extensions  parallel  to  the  three  element 
sides.  For  instance,  side  IJ 

o 


JIJ 


LIJ  ' LIJ 


‘IJ 


(15) 


where  L-^j  and  L°j  are  the  lengths  of  the  side  IJ  in  the  deformed 
and  undeformed  positions  respectively.  Since  a constant  membrane 
strain  element  is  employed  the  three  extensions  uniquely  determine 
the  membrane  strain  in  the  element. 


The  calculation  of  rotations  at  a typical  node  for  a plate 
element  involves  both  the  [E]  transformation  for  the  element  and 
the  [B]  transformation  for  the  node  and  is  based  on  the  usual  thin 
plate  assumption  that  lines  which  are  initially  straight  and  nor- 
mal to  the  surface  of  the  plate  remain  so  during  the  deformation . 
In  the  present  terminology,  an  initial  normal  to  the  plate,  e^ , 
located  at  node  I rotates  with  the  nodal  axes  (i.e.,  [x,y,z]  sys- 
tem) at  node  I and  its  components  in  the  nodal  system  remain  un- 
changed during  the  motion.  These  components  are  given  by 


In]  = [Bo]T[Eo] [e°] 


(16) 


where 

[e°]  = [0,0. 1]T 

and  [n]  are  the  components  of  the  normal  to  the  deformed  surface, 
n,  at  node  I.  It  should  be  noted  that  since  the  components  of  n 
are  constant  in  the  nodal  coordinate  system  which  rotates  during 
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the  deformation,  the  components  of  n in  the  global  or  element  sys- 
tems are  time  varying.  At  any  time  during  the  deformation  the 
angle,  9,  between  the  normal  to  the  deformed  plate  at  node  I,  n, 
and  the  current  normal  to  the  rigid  body  plane  of  the  plate,  , 
may  be  calculated  from  the  vector  cross  product 


9 


sin9 


x n 


(17) 


Assuming  that  this  is  a sufficiently  small  angle,  its  components 
in  the  element  system  are 

9 = - n (18a) 

x y v ' 

9 = n (18b) 

y x v ' 

where  the  components  of  the  deformed  normal,  n,  in  the  element 
system  are  determined  from 


[n]  = [E)T[B][n] 


(19) 


The  rotations  6 ,9 
x y 

sequent  calculations 
forces  and  moments 


(two  per  node  per  element)  are  used  in  sub- 
of  element  curvatures  and,  in  turn,  element 


2.3.2  Beam  Elements— The  element  coordinate  system  for  the 

three-dimensional  beam  element  (Figure  2)  is  constructed  somewhat 

differently  than  for  the  plate  element.  In  the  undeformed  posi- 

tion,  the  unit  vector  e^  lies  along  the  axis  of  the  beam  (the  xQ 

axis)  and  is  defined  by  the  positions  of  the  end  point  I and  J. 

The  unit  vector,  e^  (the  yQ  axis)  is  directed  normal  to  the  beam 

axis  and  lying  in  a plane  defined  by  a reference  or  pointer  node 

specified  in  the  input  data  and  containing  one  of  the  principal 

axes  of  the  beam  cross  section.  The  third  unit  vector  is  deter- 
o o o ^ 

mined  by  eQ  = e-,  x en  and  establishes  the  z axis.  In  this  case  two 
J 3 1 2 o 

unit  vectors  ri  and  £ which  initially  coincide  with  the  undeformed 

A A . 

xq  and  yQ  axes  are  embedded  in  the  nodal  coordinate  system  [x,y,z] 
at  each  end  of  the  beam.  The  nodal  components  of  these  vectors 
are 

(nl  = [Bo]T[Eo] [e°]  (20) 
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15 


where 


[e°]  = [1,0, 0]T 

and 

[C]  = [B0]T[Eq] [e°]  (21) 

where 

[e°]  = [0,1,0] 


The  construction  of  element  axes  for  the  deformed  beam  differs 
somewhat  from  the  procedure  employed  for  the  plate  since  the  po- 
sition data  for  the  two  nodes  I and  J are  not  sufficient  to  locate 
the  y and  z axes  uniquely,  although  they  are  sufficient  to  estab- 

A - \ A 

lish  the  x axes  and  the  e-^  unit  vector  directly.  The  z is  taken 
to  be  normal  to  the  plane  defined  by  the  x axis  and  the  average 
of  the  unit  vectors  and  ^ at  the  ends  of  the  beam 

e^  = e^  x (C^+r^)/2  (22) 

A 

and  the  y axis  is  formed  to  complete  the  right-handed  system 


(23) 


Displacements  in  the  element  coordinate  system  are  determined 
in  a manner  analogous  to  that  used  for  the  plate  element.  The  ex- 
tension of  the  beam  is  calculated  directly  from  the  deformed  and 
undeformed  lengths  as  determined  by  the  displacement  data  for  the 
end  points 


E 


IJ 


LIJ  ' LU 


L 


IJ 


(24) 


The  rotations  related  to  bending  at  node  I are  found  from  the  unit 
vector  n which  defines  the  deformed  position  of  the  beam  axis  at 
the  node  and  e-^  which  lies  along  the  x axis  by  means  of  the  rela- 
t ion 


9 = 


el  x 


-2k 

n 


(25) 
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(26a) 

(26b) 


For  small  angles  this  results  in 

A A 


where  the  components  of  q in  the  element  system  are  determined 
from 

In]  = [E]t[b]  [iT]  (27) 

The  torsional  rotation  0 is  determined  directly  from  the  cross 
product  relation  (^  x ZJ)  in  terms  of  components  [ n -j- ] , [ rj j. ] in 
the  element  system. 


2.3.3  Rigid  Links  —The  computer  program  incorporates  a rigid 
link  or  mas  ter- s lave  node  relationship  whereby  one  or  more  nodes  may 
be  identified  as  being  rigidly  connected  to  and  driven  by  another 
master  node.  The  slave  or  driven  nodes  have  no  mass  properties  and 
their  motions  are  determined  by  the  displacements  and  rotations  of 
the  master  node  and  the  appropriate  kinematic  relations  for  a con- 
necting rigid  link  which  is  embedded  in  the  nodal  (rigid  body)  axes 
for  the  master.  This  feature  is  intended  to  be  used  as  a means  of 
introducing  rigid  body  masses  into  a structural  model  to  simulate 
such  items  as  a vehicle  center  of  gravity,  engine  blocks,  wheel 
suspensions,  etc. 


In  terms  of  the  nodal  axes  [x,y,z]  and  transformations  [B] 
previously  defined,  the  equations  governing  these  links  can  be 
developed  in  a straightforward  manner.  If  I and  J are  taken  to 
denote  the  slave  and  master  nodes,  respectively,  then  the  position 
of  the  slave  node  with  respect  to  the  master  in  the  undeformed 
structure  may  be  written  as 


[RIG]  [ rjq]  + 


(28) 


where  [rjq]  = t^IG  ’ ^IG  ’ ZjqI  s t^ie  posi-ti-011  vector  for  node  I in 
global  coordinates  and  [A°]  is  a vector  whose  components  are  the 
initial  differences  in  position  between  the  slave  and  master  nodes. 
A rigid  body  motion  of  node  I with  respect  to  node  J implies  the 
components  of  [A^]  in  the  nodal  coordinate  system  at  node  J are 
constant  during  the  motion. 
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These  components  are 


[a]  = [bJo]t[a°]  (29) 

Thus,  at  any  time  during  the  motion  the  position  of  node  I is  de- 
termined from  the  position  of  node  J by 


[RIG1  + [uIG]  - [RJG]  +[UJG]  +[AG]  (30) 

where  [uTn]  and  [u,P]  are  the  displacements  of  nodes  I and  J and 

[A-,]  are  the  current  global  components  of  the  rigid  link  as  given 
G 

by 


[Ag]  - [Bj][A] 


(31) 


Equation  (30)  may  be  rearranged  to  give  the  displacement  of 
node  I directly  as 


^IG1  = !uJG!  + [AG]  ' [AG] 


(32) 


2.4  Element  Forces 

Having  reduced  the  displacements  to  rigid  body  and  deforma- 
tion components,  the  next  step  in  the  process  is  to  calculate  the 
forces  on  the  nodes,  taken  element  by  element,  from  the  displace- 
ments referred  to  the  element  coordinate  systems.  In  general,  the 
procedure  by  which  the  nodal  forces  are  calculated  is  derived  from 
the  principal  of  virtual  work  which  states  that  for  a kinematically 
admissible  variation  in  the  element  displacements,  the  work  done  by 
the  nodal  forces  is  equal  to  the  work  done  by  the  stresses  in  the 
element.  In  the  present  notation  this  relation  for  a typical  ele- 
ment is  written  as 

[f]T[«u]  = J [o(e)]T[6e]  dV  (33) 

V 
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wh  er  e [ f ] 

(i, j , • . • ) 


is  the  vector  of  all  forces  and  moments  at  the  nodes 
defining  the  element 


= [ f . , f . , f . ,m.  ,m.  ,m. 

L ix  iy’  iz  lx  iy’  iz 


[u]  is  the  corresponding  vector  of  element  displacements,  [a]  and 
[e]  are  the  biaxial  stresses  and  strains  at  a point  in  the  plate 


[a]  = [a 


xx 


yy  xy 


[e]  = [e 


xx 


yy 


e ] 

xy 


and  6 denotes  a variation. 


The  strains  are  related  to  the  displacements  through  a matrix 
consisting  of  the  appropriate  derivatives  of  the  element  displace- 
ment functions 


[e]  = [$(x,y ,z) ] [u] 


(34) 


Introducing  this  relation  into  equation  (33)  and  noting  that  the 
resulting  expressions  must  be  satisfied  term  by  term  in  the  virtual 
displacements,  leads  to  the  following  definition  of  the  forces 
acting  on  the  nodes. 


[ a ( e ) ] 


T[$] dV 


V 


(35) 


The  development  of  the  nodal  forces  for  particular  elements  is 
described  below  in  terms  that  are  convenient  to  each  element  formu- 
lation . 

2.4.1  Plate  Element  —The  plate  element  used  in  the  program  is 
a flat  triangular  element  lying  in  the  x-y  plane  of  the  element 
coordinate  system  defined  by  the  deformed  positions  of  the  comer 
nodes,  i,j,k  (see  Figure  3). 

/V  A /S. 

The  superscript  notation  [f],[u],  etc.  is  omitted  since  all 
quantities  involved  in  this  section  are  referred  to  a typical 
element  coordinate  system  [x,y,z]. 
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Figure  3 Triangular  Plate  Element 

The  element  formulation  is  based  on  small  deflection  linear  plate 

theory  involving  a linear  displacement  (i.e.,  constant  strain) 

field  for  the  midplane  deformations  and  a cubic  displacement  field 

for  the  plate  bending  deformations  as  given  by  Zienkiewicz  (Ref.  26). 

Node  point  degrees  of  freedom  consist  of  the  two  in-plane  motions 

u and  u which  determine  the  midplane  strains  and  the  two  nodal 
x y r 

rotations  Qx  and  9^  which  determine  the  plate  curvatures.  The  out- 
of-plane  disp lacements } uz , of  the  nodes  are  zero  since  the  x-y  ele- 
ment coordinate  plane  was  established  by  the  deformed  position  of 
the  nodes.  In  addition,  nodal  rotations  normal  to  the  plane  of 
the  plate,  0 , are  not  admitted  as  element  degrees  of  freedom  (as 
in  Ref.  26)  for  the  element  formulation  which  is  employed. 
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Strain  Displacement  Relations:  It  is  convenient  to  repre- 

sent the  element  strains  [e(x,y,z)]  in  terms  of  the  midplane 
strains  [e°(x,y)  ] and  curvatures  [ic(x,y)]  in  the  form 

[ e ] = [ e° ] - z [ k ] (36) 


where 


o 

o 

p 

o 

p 

XX  ’ 

yy  ’ 

xy 

XX  ’ 

Kyy’ 

Kxy 

T 

T 


since  these  quantities  uncouple  the  midplane  strain  and  out-of- 
plane  bending  displacement  fields . By  regrouping  the  displacements 
and  rotations  at  the  element  nodes  in  a corresponding  manner 
the  midplane  strains  and  curvatures  may  be  written  as 

[e°]  = m[u] 


(37) 


[ k ] = [ip]  [0] 


where  [$]  and  [ip]  together  are  equivalent  to  [ $ ] in  equation  (35) 
and 


[u] 

[9] 


[u.  ,u.  ,u.  ,u.  , u,  , u,  1 

L ix  ly  j x j y kx  ’ ky 

[0.  ,9.  ,9.  ,0.  ,9-|  , 9,  ] 

1 ix  ’ ly  jx  jy  kx  ky 


The  elements  of  [$]  and  [ ip ] are  derivable  from  the  shape  functions 
given  in  Ref.  26.  In  particular  [d]  is  constant  over  the  element 
and  in  terms  of  local  element  coordinates  (with  x.  =0  and  y^=0), 
has  the  form 


N] 


yk’ 

0 , 

■yr 

0 - 

0 , 

_xk  ’ 

0 , 

X. 

] 

(38) 

xk’ 

xj  ’ 

"yi  = 

where  A is  the  area  of  the  element.  The  elements  of  [ ip  ] vary  lin- 
early over  the  plate  element  and  are  obtained  by  differentiation 
of  the  cubic  shape  functions  as  given  in  Ref.  26  and  by  Meek 
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(Ref.  2 7)  with  the  curvatures  defined  as 


K 

XX 

r ,2  i 

9 uz 
9x2 

K 

yy 

92u 

z 

8y2 

K 

L xyJ 

9 u 

2 2 
_ 9x9yJ 

= m [9] 


(39) 


where  u is  the  transverse  plate  displacement.  The  derivation  of 
z 

the  elements  of  [ip]  is  given  in  detail  in  Appendix  A. 


Nodal  Forces  and  Moments : Forces  and  moments  acting  at  the 

nodes  are  calculated  from  equation  (35)  as  specialized  to  the 
nomenclature  for  this  element.  The  in-plane  forces  are  obtained 
from  the  integrals 


f.  "1 

LX 

f. 

f . 

JX 

fjy 

fkx 


kyJ 


[P]  [tf]dA 


(40) 


where  the  cross  section  plate  forces  [P]  are  obtained  by  inte- 
grating through  the  thickness 


[P]  = 
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(41) 
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L 


In  a similar  manner  the  nodal  moments  are  calculated  from  inte- 
grals of  the  form 


m . 1 

IX 
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iy 
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jx 
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jy 

“kx 
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kyj 


■l 


[MKt^dA 


(42) 


where  the  cross  section  plate  moments  [M]  are  also  obtained  by 
integration  through  the  thickness 


[M] 


rM 

XX 


h/2 

-f  Z[a(e)]dA 
-h/2 


(43) 


The  remaining  nodal  forces  f.  , f.  , f , are  found  from  overall 

3-  Z J Z KZ 

equilibrium  of  the  element  after  the  nodal  moments  have  been  es- 
tablished. Taking  moments  about  node  i,  the  origin  of  the  element 
coordinate  system,  results  in  the  following  relation  for  the 
forces 


r f . i 

JZ 

1 

'y3 

"x3 

n' 

1 , ^ 
M-l 
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2A 

:y2 

x2_ 

The  remaining  force  f^  is 

f.  = -f . - f, 

iz  jz  kz 


m.  + m . . + m, 
ly  jy  ky 

-m.  + m.  + m,  , 

ix  j x kx  J 


(44) 


found  from  transverse  equilibrium 

(45) 


Numerical  Integrations : The  integrals  involved  in  the  defi 

nitions  of  cross  section  forces  and  moments  (equations  (41)  and 
(43))  are  evaluated  numerically  at  the  midpoints  of  the  three 
sides  of  the  triangle  (points  1,  2,  and  3 in  Figure  3).  The 
integration  procedure  involves  a piecewise  linear  trapezoidal 
integration  rule  employing  stresses  at  'n'  evenly  spaced  points 
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through  the  thickness  with  the  first  and  last  points  lying  on  the 
outer  surfaces  of  the  plate.  At  least  two  points  must  be  used 
since  this  number  provides  exact  results  for  elastic  problems  and 
four  or  five  points  have  proven  to  provide  satisfactory  accuracy 
for  all  cases  in  which  detailed  accuracy  studies  have  been  made. 
The  general  forms  given  in  equations  (41)  and  (43)  apply  to  any 
consistent,  nonlinear  biaxial  stress  strain  relations  although 
applications  to  date  have  been  limited  to  bilinear  elastic-plastic 
material  properties.  It  should  be  noted  that  the  formal  defini- 
tion of  internal  forces  and  moments  has  been  maintained  in  the 
computer  program  so  that  explicit  relations  among  extension,  cur- 
vature, force  and  moment  may  be  incorporated  with  relative  ease. 
This  has  been  found  useful  in  purely  elastic  problems  as  a means 
of  improving  the  efficiency  of  the  program. 

Integrations  over  the  surface  of  the  plate  (equations  (40) 
and  (42))  are  evaluated  by  three  point  Gaussian  integration  using 
the  three  midside  points  described  above.  In  this  case  (see 
Ref.  26)  the  weighting  factors  are  1/3  and  a typical  integration 
(equation  (42)  for  instance)  results  in 


where  [M^]  and  [i|^]  are  the  values  at  the  midside  points. 

Stress -Strain  Relations:  The  stress- strain  relations  em- 

ployed in  the  plate  element  formulation  are  those  of  a bilinear, 
elastic-plastic  strain  hardening  material  as  indicated  in  Figure  4. 
The  plastic  flow  calculations  are  based  on  a Mises  yield  criterion 
and  isotropic  hardening  for  the  plane  stress  conditions  which  are 
assumed  to  apply  in  the  plates.  The  computational  algorithm  fol- 
lows that  presented  by  Hartzman  and  Hutchinson  (Ref.  28)  as  spe- 
cialized for  small  strain  and  plane  stress  conditions.  Given  a 


small  increment  in  strain  from  the  prior  to  the  current  state 

Ae  , Ae  , Ae  the  current  stress  state  is  established  by  the 
xx  yy  xy  J 

following  procedure. 


A 


L = 1,2,3 


(46) 


prior  stress  state  at  a point  in  the  plate,  a 
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Figure  4 Uniaxial  Bilinear  Stress-Strain  Relation 


A tentative,  current  stress  state 
though  the  strain  increment  were 


axx’ ayy ’ axy  is  calculated  as 
completely  elastic: 
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(47) 

(47) 


where  E,  v and  G are  the  usual  elastic  constants.  The  correspond 
ing  value  of  effective  stress  for  the  Mises  yield  criterion  is 
also  calculated  from  the  expression 
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- - -2  3 - 

a a + a + T a 

xx  yy  yy  2 xy 


1/2 


(48) 
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which  is  Used  to  determine  if  plastic  flow  has  taken  place  during 
the  strain  increment.  If  oe  is  less  than  the  most  recent  effec- 
tive stress  q°  at  which  yielding  occurred  the  material  is  unload- 
ing and  the  tentative  stresses  are,  in  fact,  the  current  stress 

values  a a ,a  . Otherwise  the  tentative  stress  state  must  be 
xx’  yy  xy  . A D ^ OQ 

modified  to  correct  for  plastic  behavior.  As  shown  m Ret.  Z8 , 

the  true  value  of  the  current  effective  stress  will  be 
[°e  + (3G)  °e ] 


°e 


1 + (Tg) 


(49) 


where  H'  is  the  slope  of  the  stress  versus  plastic  strain  curve 
given  by 

E E 

H.  = E. 

H E + E 

P 

with  Ep  being  the  plastic  modulus  shown  in  Figure  4.  The  current 
stress  state  is  then  given  by 


a = 7i — \ ta  + X (a  +a  „)] 

xx  (1  + 3X)  xx  xx  yy 

a = TT — [a  + X(a  +0  )] 

yy  (1  + 3X)  yy  xx  yy 

1 r-  i 

o -i  I q L o J 

xy  1 + 3X  xyJ 

where  the  factor,  X,  is  given  by 


(50) 


Further  details  concerning  the  development  of  this  algorithm  may 
be  found  in  Ref.  28. 


2.4.2  Beam  Element— The  beam  element  employed  in  the  program 
is  based  on  the  conventional  small  deflection  formulation  involv- 
ing cubic  displacement  fields  (i.e. , linear  curvature  distribution) 
for  the  transverse  displacements  and  linear  displacement  distribu- 
tions for  the  axial  and  torsional  motions  (i.e.,  constant  axial 
strain  and  twist).  The  development  of  the  element  forces  paral- 
lels that  used  for  the  plate  with  such  changes  and  simplifications 
as  are  appropriate  for  beam  theory. 
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Strain -Displacement  Relations:  The  longitudinal  strain  in 

the  beam,  e „ (x) , is  written  in  terms  of  the  extension  at  the 

XX  ' 

centroid  e°(x)  and  curvatures  k (x)  and  k (x)  about  the  principal 

yy  zz 

axes  of  the  cross  section  (see  Figure  5)  in  the  form 
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yy 


(51) 


The  extensions  and  curvatures  are  in  turn  written  in  terms  of  the 
nodal  displacements  using  the  appropriate  derivatives  of  the  dis- 
placement fields  which  results  in 
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(52) 


where  L is  the  length  of  the  beam  and  £ = X/L. 

Nodal  Forces  and  Moments:  Forces  and  moments  acting  at  the 
nodes  are  determined  by  introducing  these  strain  displacement  re- 
lations, equations  (51)  and  (52),  into  the  appropriate  form  of 
equation  (35)  with  the  following  results: 
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Figure  5 Seam  Coordinates  and  Cross  Section  Geometry 
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where  M is  the  internal  cross  section  moment  about  the  y axis 

yy 

M = / z a dx 

yy  J xx 

A 

“zi  - lM'2)  Mzz  <55> 

L 

m = r I (H-l)  M dx 
zy  LJ  zz 

L 

where  M is  the  internal  cross  section  moment  about  the  z axis 
z z 

M = - y a dA 
zz  J J xx 

A 

The  torsional  moments  are  omitted  from  the  formal  development 
given  above  and  are  obtained  directly  in  terms  of  the  tor- 
sional stiffness  for  the  beam  section  as 
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where  is  the  torsional  constant  for  the  cross  section  and  G is 
the  shear  modulus.  Finally,  the  transverse  nodal  forces  are  de- 
termined from  equilibrium  conditions  using  the  nodal  moments  as 


m . + m 

JL 1 Z1 
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m . +m  . 
zi  zj  . 


(57) 


Numerical  Integrations:  The  integrals  involved  in  the  defi- 

nition of  internal  cross  section  moments  and  forces  are  obtained 
by  numerical  integration  of  the  stress  distribution  over  the  cross 
section.  The  cross  sections  are  treated  as  thin-walled  plate  work 
sections  defined  as  indicated  in  Figure  5 by  tables  of  the  coor- 
dinates of  points  on  the  centerline  of  the  wall  and  the  thicknesses 


The  definitions  of  the  cross  section  moments,  Myy  and  Mzz  are 
chosen  to  be  in  the  positive  sense  on  the  positive  face  of _ the 
beam  cross  section  and  to  result  in  the  usual  sign  convention  for 
elastic  material,  i.e.,  Mzz = E !zz  Kzz  and  Myy  = ~EIyy  KYY' 
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of  wall  sections.  The  integration  algorithm  is  a generalization 
to  two  dimensions  (the  y-z  plane)  of  the  procedure  employed  in  the 
plate  element  and  assumes  a piecewise  linear  variation  in  stress 
around  the  periphery  of  the  cross  section.  As  in  the  plate  cross 
section  integrations  this  procedure  preserves  the  exact  elastic 
stiffness  of  the  cross  section  provided  that  integration  points 
include  all  corners  in  the  cross  section.  The  stress  strain 
algorithm  relation  employed  in  these  integrations  is  the  one- 
dimensional version  of  that  prescribed  in  connection  with  the 
plate  element  which  can  be  obtained  from  equations  (47)  through 
(50)  simply  by  deleting  all  occurrences  of  c ^ and  aXy-  As  with 
the  plate  the  definitions  of  internal  cross  section  forces  and 
moments  have  been  maintained  in  the  computer  program  to  facilitate 
subsequent  use  of  explicit  definitions  of  these  quantities. 

The  integrals  over  the  length  of  the  beam  are  not  carried 

out  explicitly  but  instead  are  evaluated  in  terms  of  P , M , M 

xx  yy 

at  the  end  points  of  the  beam: 
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(58) 


For  the  moments  this  is  equivalent  to  assuming  that  the  internal 
moments  M and  M zz  vary  linearly  along  the  beam  and  is  consistent 
with  the  usual  beam  formulation  which  implies  constant  shear  and 
linear  moment  variation  along  the  beam 


2.5  Nodal  Force  Transformations 

The  preceding  section  describes  the  manner  in  which  forces 
and  moments  acting  on  the  nodes  of  a given  element  are  calculated 
for  specified  values  of  nodal  displacement.  These  are  calculated 
on  an  element  by  element  basis  with  reference  to  the  individual 
element  coordinate  system  (x,y,z)  and  must  be  transformed  to  the 
global  (Xg.y^.Zg)  and  nodal  (x,y,z)  systems  respectively  in  order 
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to  be  entered  into  the  appropriate  equations  of  motion  (equations 

/\ 

(1)  and  (2)).  For  a particular  element  the  nodal  forces,  [f], 
are  transformed  to  the  global  system  by  means  of  the  element 
transformation 

[fG]  = [E] [f]  (59) 

and  entered  directly  into  the  force  summations  in  equation  (1). 

The  moments  [m]  are  first  transformed  to  the  global  system  and, 
if  the  nodes  in  question  are  not  part  of  a master-slave  relation, 
subsequently  to  the  nodal  coordinate  system  for  the  particular 
node 

[m]  = [ B ]T  [ E ] [m]  (60) 

and  entered  into  the  moment  summations  (equations  (2)).  If  the 
node  in  question  is  associated  with  a master  node,  the  moments 
are  first  transformed  to  the  global  system  by  the  transformation 

[mGI]  = [E] [mx]  (61) 

where  I denotes  quantities  at  the  slave  node,  and  then  transferred 
to  the  master  node  (denoted  by  J)  by  the  cross-product  relation 

mJG  = mIG  + AG  X fG 

The  moments  for  the  master  node  [m^p]  are  then  transformed  to  the 
nodal  system  for  that  node 

[mj]  = [B]T[mJG]  (63) 

and  entered  into  the  appropriate  moment  summations  in  equation  (2) . 
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3.  ANALYSIS  RESULTS 


A substantial  number  of  test  and  demonstration  problems  have 
been  treated  during  the  development  of  the  analysis  and  computer 
program  reported  here.  For  the  most  part  the  test  problems  dealt 
with  relatively  simple  structural  configurations  having  known  re-t 
suits  and  were  selected  to  isolate  and  test  particular  aspects  of 
the  formulation  and  program  during  development.  Thirteen  such 
test  problems  of  successively  increasing  complexity  were  formally 
reported  in  progress  reports . In  all  cases  the  analyses  were  pur- 
sued until  satisfactory  comparison  with  previous  results  were  ob- 
tained or  at  least  to  the  point  at  which  differences  between  re- 
sults could  clearly  be  reconciled  on  the  grounds  of  differences  in 
loadings,  boundary  conditions,  or  ambiguities  in  the  results  used 
for  comparison  purposes.  The  contents  of  these  test  problems  are 
recapitulated  briefly  in  Table  1,  and  results  for  two  more  com- 
plicated problems,  namely  dynamic  snap-through  for  an  elastic 
dome  and  the  response  of  an  elastic-plastic  panel  under  impulsive 
load  are  reviewed  in  the  following  sections.  Results  are  also 
presented  in  subsequent  sections  which  involve  the  simulation  of 
two  actual  crash  events  taken  from  a previous  series  of  full  scale 
barrier  crash  tests. 


3.1  Test  Problems 

Results  for  two  test  problems  are  presented  in  this  section. 
These  have  been  chosen  to  demonstrate  the  capabilities  of  the  an- 
alysis procedure  in  terms  of  fairly  difficult  problems  for  which 
some  results  are  available  in  the  literature.  It  should  be  noted 
however  that  the  results  used  as  references  are  also  based  on  nu- 
merical computations  and  approximations  and  are  not  considered  to 
be  "exact"  in  drawing  comparisons. 

3.1.1  Dynamic  Buckling  of  a Spherical  Shell  Segment  — This 
problem  concerns  the  dynamic  response  of  a clamped,  axisymmetric , 
elastic  spherical  cap  under  suddenly  applied  external  pressure. 
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SUMMARY  OF  TEST  PROBLEMS 
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As  indicated  in  Figure  6,  the  cap  consists  of  a segment  of  a sphere 
of  60  in.  radius,  with  a shallow  rise  of  0.30  in.  and  a thickness 
of  0.08  in.  A plan  view  of  the  finite  element  mesh  employed  in 
the  solution  is  shown  in  Figure  7.  The  mesh  consists  of  144  plate 
elements  and  82  nodes  and  covers  the  complete  dome  with  no  con- 
straints or  geometric  features  introduced  to  enforce  an  axisymmetric 
response  pattern.  This  problem  has  also  been  treated  by  Archer 
and  Lange  in  Ref.  30  by  means  of  a numerical  integration  of  the 
finite  difference  equations  derived  from  Marguerre's  equations  for 
the  axisymmetric  response  of  spherical  shells.  Displacement  time 
histories  for  the  vertical  motion  of  the  center  of  the  cap  under  an 
inward  pressure  of  44  psi  are  shown  in  Figure  6 for  the  present  re- 
sults and  for  Ref.  30.  It  is  evident  that  there  is  fairly  close 
agreement  between  the  two  results  with  the  finite  element  response 
being  slightly  delayed  (i.e.,  stiffer)  than  the  results  of  Ref.  30. 
It  is  also  apparent  that  the  response  lies  in  the  snap-through 
range  since  the  maximum  displacement  approaches  twice  the  unde- 
formed rise  of  the  shell.  This  is  also  clearly  shown  in  the  suc- 
cessive deformed  shapes  of  the  shell  shown  in  Figure  8. 

3.1.2  Dynamic  Response  of  a Cylindrical  Panel— The  purpose  of 
this  test  problem  is  to  compare  the  results  of  the  present  analy- 
sis against  previously  obtained  results  involving  dynamic,  elastic- 
plastic  shell  buckling  response.  As  indicated  in  Figure  9,  the 
problem  consists  of  a 120  deg  segment  of  an  aluminum  cylindrical 
shell,  supported  and  clamped  at  all  four  edges  and  loaded  explosive- 
ly over  a portion  of  its  surface,  as  shown.  Previous  experimental 
and  analytical  results  are  available  for  this  problem  in  the  liter- 
ature, particularly  in  Refs.  31  and  32. 

The  mes,h  geometry  employed  is  illustrated  in  Figure  10  and 
consists  of  64  elements  and  45  nodes  representing  one-half  of  the 
symmetrical  shell.  The  disp lacement- time  history  of  node  5 is 
shown  in  Figure  11  together  with  corresponding  results  from  Refs. 

31  and  32.  The  explosive  load  according  to  Refs . 31  and  32  corre- 
sponds to  an  initial  (i.e.,  impulsive)  velocity  of  5650  in. /sec 
directed  radially  inward  over  the  loaded  portion  of  the  slid  11 . 
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Figure  6 Displacement-Time  History  of  a Clamped  Spherical  Shell 
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Figure  7 Spherical  Shell  Finite  Element  Mesh 
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Explosive  Load 


a.  Geometry 


° AL  6Q61-T6 

E = 10.5  x 106  psi;  v = 0.333 
p = 2.54x10  ^ lbf  - sec^/in 

o = 44,000  psi;  E =0 
y - p 

b . Properties 

injure  9 Cylindrical  Shell  - Geometry  and  Properties 


Figure  10  Cylindrical  Shell  Segment  - Undeformed 


Displacement  (in. 


Time  (msec) 

Figure  11  Cylindrical  Shell  Segment  Crown 

Displacement  versus  Time  (Node  6) 
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This  corresponds  to  the  constant  velocity  line  in  Figure  11  which 
is  tangent  to  the  displacement  curve  at  t = 0.  The  present  re- 
sults are  bracketed  by  those  of  Refs.  31  and  32  with  the  improved 
analysis  results  of  Ref.  32  lying  below  and  the  experimental  re- 
sults above.  The  older  analysis  results  from  Ref.  31  are  very 
close  to  the  experimental  curve  but  subsequent  improvements  in 
that  procedure  (Ref.  32)  result  in  lower  displacements.  Improved 
correlations  were  obtained  in  Ref.  32  by  introducing  yielding 
boundaries  and  this  was  not  attempted  in  the  present  work. 

The  first  part  of  the  velocity  history  for  node  5 is  shown 
in  Figure  12.  The  velocity  decreases  from  its  initial  value  of 
5650  in/sec , increases  again  as  the  shell  rebounds,  and  then  de- 
creases fairly  smoothly  to  zero  at  approximately  0.4  msec.  The 
corresponding  displacement  pattern  at  0.4  msec  is  shown  in  Fig- 
ure 13.  While  further  efforts  could  be  made  to  improve  the  cor- 
relations shown  in  Figure  11,  principally  in  the  boundary  condi- 
tions and  in  a refined  mesh,  the  present  results  are  considered 
satisfactory . 

3.2  Barrier  Test  Simulations 

The  test  problems  presented  throughout  the  progress  reports 
and  in  the  preceding  section  are  an  essential  means  of  assuring 
that  the  analysis  formulation  and  computer  program  are  capable  of 
describing  the  types  of  phenomena  for  which  they  are  intended. 
However,  the  transition  from  this  class  of  problems,  which  are 
fairly  small  and  idealized,  to  the  analysis  of  actual  vehicle 
structures  in  realistic  crash  events  is  by  no  means  a trivial  step. 

In  the  latter  case  the  analysis  must  deal  with  much  longer  dura- 
tion events  and  structures  of  substantially  greater  size  and 
geometric  complexity  with  realistic,  rather  than  idealized,  bound- 
ary conditions  and  connection  details.  In  fact,  the  proper  and 
efficient  application  of  an  analysis  procedure  to  real  structures 
and  events  is  nearly  as  much  a research  undertaking,  at  least 
initially,  as  is  the  original  formulation  and  coding  of  the  analysis. 
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Velocity  (in. /sec) 
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Figure  13  Cylindrical  Shell  Segment  - Deformed 


Also,  the  planning  and  execution  of  such  analyses  has  many  as- 
pects in  common  with  the  planning  and  execution  of  full  scale 
crash  tests,  not  the  least  of  which  is  the  acquisition  and  hand- 
ling of  the  large  quantities  of  data  which  are  inherent  in  either 
activity. 

This  section  describes  the  results  of  our  efforts  to  obtain 
relatively  large  scale  simulations  of  actual  vehicle  crash  events. 
The  events  which  are  the  objects  of  these  simulations  were  se- 
lected from  a series  of  barrier  tests  conducted  by  Dynamic  Science 
(Ref.  33).  The  series  of  tests  of  interest  involved  frontal  im- 
pacts at  various  velocities  of  a 1968  Plymouth  Fury  having  differ- 
ent portions  of  the  front  end  structure  and  machinery  removed  in 
different  tests  as  summarized  in  Table  2.  The  object  of  these 
tests  was  to  isolate  and  identify  the  dynamic  resistance  proper- 
ties of  the  various  front  end  components  in  three  groups,  namely: 
frame  structure,  sheet  metal  structure  and  the  engine  and  machinery 
components.  Simulations  were  developed  for  two  test  configurations 
involving  the  stub  frame  alone  (test  283-51)  and  the  stub  frame 
and  sheet  metal  together  (test  283-53)  since  these  configurations 
did  not  involve  the  added  complications  associated  with  the  engine 
components.  The  test  involving  sheet  metal  alone  was  excluded  due 
to  the  extreme  and  unrealistic  deformation  patterns  which  appeared 
in  the  test  results. 


Table  2 

PLYMOUTH  FURY  TEST  SERIES 


Test 

No. 

Configuration 

Total 

Weight 

(lbs) 

Imp  ac  t 
Veloci ty 
(mph) 

283-49 

Sheet  metal 

2665 

19.85 

283-50 

Engine 

3914 

19. 48 

283-51 

Stub  frame 

2909 

29.64 

283-52 

Full  vehicle 

4192 

51.45 

283-53 

Stub  frame  and  sheet  metal 

2997 

43.99 
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The  finite  element  models  developed  for  these  two  test  con- 
figurations represent  one-half  of  the  structural  components  for- 
ward of  the  vehicle  fire  wall  and  are  based  on  an  assumption  of 
symmetry  for  both  the  vehicle  and  its  response  modes  about  the 
center  plane.  The  frame  model  consists  entirely  of  beam  elements 
and  the  frame  and  sheet  metal  model  contains  the  frame  model,  plate 
elements  representing  the  sheet  metal,  and  additional  beam  elements 
representing  stiffener  sections  in  the  hood.  In  both  cases  the 
appropriate  inertial  properties  of  the  remainder  of  the  vehicle 
are  provided  as  added  masses  at  the  vehicle  center  of  gravity  and 
are  attached  to  the  frontal  structure  by  means  of  rigid  links. 

Data  relating  to  the  tests  and  the  vehicle  geometry  were  obtained 
from  several  sources  including  the  following: 

• Reports  and  high-speed  motion  pictures  associated  with 
the  test  program  (obtained  through  DOT) 

• Sheet  metal  profiles,  thicknesses,  and  material  prop- 
erties obtained  from  Stanford  Research  Institute 

• Additional  measurements  and  observations  of  a 1968 
Plymouth  Fury  by  IITRI  staff. 

3.2.1  Stub  Frame  Impact  (Test  283-51)  — Plan  and  side  views  of 
the  finite  element  mesh  for  the  stub  frame  are  shown  by  the  dashed 
(i.e.,  undeformed)  lines  in  Figures  14  and  15.  The  model  consisted 
of  22  nodes  and  21  beam  elements  plus  a rigid  link  connecting  the 
rear  of  the  frame  at  the  firewall  attachment  point  to  the  center 
of  gravity  of  the  vehicle.  Additional  mass  was  lumped  at  the  ap- 
proximate location  of  the  tire  connection  and  also  at  the  vehicle 
center  of  gravity  in  order  to  introduce  the  proper  inertial  charac- 
teristics into  the  analysis.  Particular  emphasis  was  placed  in 
accurately  modeling  the  various  open  and  closed  wall  beam  sections 
that  comprise  the  stub  frame.  Six  different  beam  cross  sections 
were  employed  in  the  model  as  summarized  in  Table  3.  In  all  cases 
the  material  properties  were  taken  as  those  of  mild  steel  with  an 
elastic  modulus  of  30x10  psi,  a Poisson's  ratio  of  0.3  and  a 
yield  stress  of  36,000  psi.  Boundary  conditions  consisted  of 
fixing  the  fore  and  aft  translatory  motion  of  the  nodes  in  imme- 
diate contact  with  the  barrier. 
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Direction  of  Impact 
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Figure  15  30  mph  Stub  Frame  Impact  at  40  msec  (Side) 
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Rebounding  of  the  frame  from  the  rigid  barrier  was  not  taken 
into  account,  therefore,  once  a node  came  in  contact  witu  the 
barrier,  it  remained  in  contact.  It  should  be  noted  that  at  the 
beginning  of  the  solution  phase  of  the  analysis  it  was  assumed 
that  all  nodes  along  the  bumper  were  in  contact  with  the  barrier. 

The  structure  was  driven  with  an  initial  velocity  of  30  mph  at 
all  nodes. 

Top  and  side  views  of  the  deformed  frame  at  40  msec  after 
impact  are  also  shown  as  the  solid  lines  in  Figures  14  and  15.  It 
should  be  noted  that  the  deformed  structures  are  plotted  to  the 
same  scales  as  the  undeformed  mesh  and  no  magnification  factors 
have  been  introduced.  In  the  plan  view,  Figure  14.  all  members 
show  substantial  forward  motion  with  bowing  in  the  forward  longi- 
tudinal frame  section  and  collapse  of  the  bumper  supports.  Simi- 
lar behavior  is  seen  in  the  test  films  including  bowing  of  the  frame 
members  and  severe  distortions  in  the  bumper  supports.  In  the 
films  the  triangle  formed  by  the  bumper  supports  also  pivots  away 
from  the  wall  after  collapse  but  this  was  prevented  by  the  bound- 
ary constraints  used  in  the  analysis.  In  the  side  view,  Figure  15, 
the  frame  bends  throughout  the  forward  portion  containing  the 
bumper,  bumper  supports,  and  forward  frame  are  seen  to  have  climbed 
the  barrier  a distance  of  approximately  15  in.  This  same  pattern, 
and  particularly  the  upward  motion  of  the  front  part  of  the  frame, 
is  clearly  evident  in  the  test  films  although  in  the  films  the 
collapse  of  the  bumper  supports  is  somewhat  more  severe  and  the 
end  of  the  forward  frame  section  comes  in  contact  with  the  wall. 

At  this  point  in  the  motion  the  frame  has  begun  to  rebound  from 
its  maximum  bent  position  and  the  center  of  mass  begins  to  pivot 
downward  about  the  point  of  contact  with  the  wall.  At  this  point 
the  film  sequences  indicate  that  this  motion  is  reacted  by  the 
rear  wheel  suspension  systems  but  in  the  analysis  this  pivoting 
continues  since  the  rear  wheels  were  not  incorporated  in  the  model, 
nor  were  any  vertical  constraints  applied  to  the  center  of  gravity. 
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Longitudinal  velocity  and  acceleration  at  the  vehicle  center 
of  gravity  for  40  msec  of  motion  from  the  analysis  and  test  re- 
sults are  shown  in  Figures  16  and  17.  As  seen  in  the  velocity- 
time histories,  the  vehicle  has  lost  only  about  20  percent  of 
its  initial  impact  velocity  (about  36  percent  of  its  initial  ki- 
netic energy).  In  this  test  however,  this  time  period  contains 
the  maximum  vehicle  decelerations  and  barrier  loads.  A time 
shift  of  approximately  10  msec  between  the  analysis  and  test  re- 
sults is  clearly  evident  in  the  velocity  histories  which  may  be 
the  result  of  filtering  in  the  test  data  or  differences  in  zero 
time  between  the  test  and  analysis  procedures.  Whatever  the 
cause,  this  time  shift  is  found  in  comparing  all  corresponding 
results  and  will  not  be  commented  on  further  in  subsequent  discus- 
sion. The  acceleration  time  record  obtained  from  the  analysis 
(Figure  17)  shows  an  early  deceleration  peak  of  about  10  g in  the 
first  3 msec  (note  that  a signal  transit  from  bumper  to  center  of 
gravity  takes  approximately  0.25  msec)  and  a series  of  sharp  peaks 
on  the  order  of  30  g between  15  to  27  msec.  Peaks  in  the  data  are 
about  22  g at  12  msec  and  40  g at  30  msec  with  a substantial  re- 
bound peak  at  20  msec.  In  this  case  the  test  data  clearly  show 
the  effects  of  filtering.  The  early  peak  in  the  analysis  clearly 
is  related  to  the  initial  resistance  and  collapse  of  the  bumper 
supports  while  the  subsequent  peaks  are  probably  attributable  to 
collapse  and  overall  bending  of  the  frame.  The  analysis  also  rou- 
tinely provides  contact  forces  at  the  barrier  but  in  this  case 
these  data  were  lost  beyond  the  first  5 msec  of  resoonse. 

As  indicated  in  Figure  18  however  the  early  value  of  the  total 
barrier  force  was  approximately  23,000  lbs  (for  half  of  the  ve- 
hicle) as  compared  to  a value  of  about  45,000  lbs  in  the  corre- 
sponding data  (for  the  full  vehicle) . This  apparent  agreement 
may  not  be  precise  however  when  filtering  in  the  data  is  taken 
into  account. 

3.2.2  Stub  Frame  and  Sheet  Metal  Impact  (Test  283-53)  - An  i so- 
me trie  view  of  the  finite  element  mesh  for  the  stub  frame  and  sheet 
metal  model  before  impact  is  shown  in  Figure  19. 
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Figure  16  Vehicle  Center  of  Gravity  Velocity  for  30  mph  Stub  Frame  Test 


Time  (msec) 


Figure  18 


Stub  Frame  Test  - Barrier  Force 
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Figure  19  Stub  Frame  and  Sheet  Metal  Model  (Undeformed) 
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Figure  20  Stub  Frame  and  Sheet  Metal  Model  (Deformed,  20  msec) 


This  model  comprises  the  left  front  quarter  of  the  1968  Plymouth 
employed  in  the  test.  The  major  components  of  the  model  are  the 
hood,  left  front  fender,  wheel  well  and  frame  and  radiator  support. 
The  wheel  well  is  shown  separately  for  the  sake  of  clarity  in  the 
figure  but  is  actually  an  integral  part  of  the  model  and  is  at- 
tached along  the  lines  defined  by  the  match  points  1 through  4 in- 
dicated in  the  figure.  The  thickness  of  all  sheet  metal  elements 
was  taken  as  0.039  in.  material.  The  frame  portion  of  the  model 
is  identical  to  that  employed  in  the  stub  frame  simulation  with  the 
addition  of  the  radiator  support  structure  indicated.  The  hood  por- 
tion of  the  model  is  free  to  move  relative  to  the  fender  at  their 
line  of  intersection,  and,  in  addition  to  the  plate  elements  shown, 
contains  2.0x1.0x0.05  in.  hat  section  stiffeners  at  reinforce- 
ments. As  with  the  stub  frame  model  additional  mass  was  lumped  at 
the  approximate  lo, cation  of  the  wheel  suspension  and  also  at  the 
vehicle  center  of  gravity.  In  this  case  all  points  along  the  rear 
edge  of  the  fender  and  hood  as  well  as  the  end  of  the  frame  were 
attached  to  the  vehicle  center  of  gravity  by  means  of  rigid  links. 
The  model  was  driven  by  an  initial  longitudinal  velocity  of  774  in./ 
sec  for  all  nodes  with  longitudinal  constraints  applied  at  nodes  in 
the  forward  section  of  the  model  successively  as  they  came  into 
contact  with  the  barrier. 

The  deformed  shape  of  the  model  at  20  msec  after  impact  is 
snown  in  Figure  20.  At  this  point  the  forward  motion  of  the  center 
of  gravity  is  approximately  15  in.  and  the  forward  portions  of  the 
frame  assembly  and  the  first  two  lines  of  nodes  in  the  hood  and 
fender  are  in  contact  with  the  barrier.  Severe  distortions  are  in- 
dicated in  the  forward  portions  of  the  model  with  the  bumper  and 
bumper  brackets  collapsed  downward  and  the  hood  and  fender  collaps- 
ing upward  and  outward  respectively.  Separation  at  the  hood  and 
fender  junction  line  is  cleaily  evident  and  the  rear  portion  of 
the  hood  has  begun  to  buckle  upward  away  from  the  fender.  The 
wheel  well  is  not  greatly  involved  although  the  forward  wall  has 
bent  to  conform  to  the  distortions  in  the  forward  section  of  the 
fender . 
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A comparison  of  longitudinal  acceleration  time  histories  at 
the  vehicle  center  of  gravity  is  given  in  Figure  21  for  zO  msec 
of  motion.  The  analysis  result  is  similar  to  that  obtained  from 
the  analysis  of  the  stub  frame  configuration  in  that  an  early 
acceleration  pulse  of  approximately  10  g associated  with  collapse 
of  the  bumper  support  system  is  apparent.  Further  the  maximum  ac- 
celeration levels  during  this  phase  of  the  response  is  approxi- 
mately 30  g as  in  the  results  for  the  stub  frame.  In  this  case 
however  the  30  g response  is  clearly  associated  with  the  impact 
of  fender  and  hood  sheet  metal  rather  than  with  substantial  frame 
bending.  The  corresponding  test  results  indicate  a peak  of  about 
25  g at  a somewhat  later  time  and  also  clearly  show  the  effects 
of  filtering  in  comparison  to  the  analysis  data. 

A similar  comparison  of  total  barrier  force  is  provided  in 
Figure  22  in  which  50  percent  of  the  test  total  is  compared  with 
the  current  result  obtained  from  a model  of  half  of  the  vehicle. 

The  analysis  results  show  an  early  force  level  of  approximately 
40,000  lbs  associated  with  the  collapse  of  the  bumper  support  and 
frame  system  which  gradually  decays  for  the  remainder  of  the  re- 
sponse. The  impact  of  the  forward  edges  of  the  hood  and  fender  ap- 
pear as  a sudden  rise  in  barrier  force  at  about  6 msec  to  a level 
of  approximately  80,000  lbs  which  then  decays  over  the  following 
10  to  15  msec.  This  peak  is  associated  with  only  the  first  line  of 
nodes  in  the  hood  and  fender  since  forces  from  the  subsequent  im- 
pact of  succeeding  sheet  metal  nodes  were  lost  due  to  a data  process- 
ing error.  In  fact,  the  first  three  lines  of  nodes  in  the  fender 
and  hood  were  involved  with  the  barrier  and  are  included  in  the  re- 
sults given  in  Figure  21.  In  Figure  22,  the  impact  of  succeeding 
nodes  would  appear  as  subsequent  peaks  in  the  force  time  history 
which  would  appear  about  the  time  that  the  peak  due  to  the  previous 
nodes  had  decayed  substantially.  In  reality,  the  progressive  con- 
tact of  the  sheet  metal  with  the  barrier  is  a more  or  less  continu- 
ous process  and  the  successive  impacts  of  the  nodes  would  amount 
to  a stepwise  approximation  of  this  process . 
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Figure  21  Stub  Frame  and  Sheet  Metal  Vehicle  Center  of  Gravity  Acceleration 
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A study  of  the  force- time  histories  for  individual  nodes  at  the 
front  of  the  hood  and  fender  indicates  that  of  the  80,000  lbs 
total  force  in  the  first  peak,  about  46,000  lbs  is  derived  from 
collapse  of  the  hood  and  about  34,000  lbs  from  the  fender. 

The  test  simulations  presented  here  are  by  no  means  a com- 
plete or  exhaustive  study  even  for  the  limited  series  of  tests 
which  were  used  for  comparison  purposes.  For  a variety  of  reasons, 
the  results  given  apply  only  to  the  initial  phases  of  the  respec- 
tive crash  events  and  improvements  could  also  be  made  in  the  models 
(the  provision  of  rear  wheels  or  other  vertical  constraints  at  the 
center  of  gravity,  for  instance)  and  in  the  associated  data  process- 
ing. Despite  these  limitations,  however,  the  results  appear  to  be 
a representative  portrayal  of  the  phenomena  associated  with  vehicle 
barrier  collisions  generally,  and  a reasonable  correlation  with  the 
particular  test  data  for  which  comparisons  were  attempted. 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  report  presents  the  results  of  a research  project  under- 
taken by  IITRI  for  DOT  to  develop  a finite  element  computer  program 
for  use  in  the  dynamic  analysis  of  vehicle  structure,  including 
sheet  metal,  in  a crash  environment.  This  research  program  involved 
the  following  major  tasks: 

• A technique  was  developed  for  the  finite  element  analy- 
sis of  the  dynamic  response  of  plate  beam  structures 
involving  very  large  displacements  and  rotations  and 
elastic -plastic  material  behavior.  The  principal  fea- 
ture of  this  technique  involves  the  decomposition  of 
the  element  displacement  field  into  rigid  body  compo- 
nents and  deformation  components  thus  allowing  the  use 

of  a small  deflection  element  formulation  in  the  analysis. 

• A computer  program  was  developed  incorporating  this 
analysis  procedure  for  beam  and  plate  elements  and  rigid 
links  together  with  appropriate  time  integration  pro- 
cedures and  material  property  descriptions. 

• A substantial  number  of  test  and  demonstration  problems 
were  analyzed  with  this  computer  program  ranging  from 
simple  classical  solutions  for  beams  and  plates  through 
large  scale  simulations  of  actual  crash  tests . 


4.1  Conclusions 

Based  on  the  work  carried  out  in  this  research  program  our 
principal  conclusions  are  as  follows: 

(1)  The  analysis  formulation,  while  somewhat  novel,  is  a for- 
mally sound  procedure  and  a proper  and  useful  approximation  technique 
for  the  dynamic  analysis  of  beam  and  plate  structures  involving  large 
deflections  and  rotations. 

(2)  The  computer  program  developed  is  a correct  rendering  of 
this  analysis  technique  and  provided  accurate  results  with  reason- 
able efficiency  in  comparison  to  available  known  solutions. 

(3)  Finally,  and  most  importantly,  the  computer  program  is  read- 
ily applicable  to  realistic  vehicle  structures  and  provides  credible 
results  for  actual  crash  events  as  demonstrated  by  the  simulations 

of  the  end-on  barrier  tests. 
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4.  2 Recommendations 


In  our  opinion  the  analysis  and  computer  program  resulting 
from  this  research  is  a substantial  improvement  over  previous 
attempts  at  a finite  element  analysis  of  vehicle  structures  and 
is  a potentially  valuable  addition  to  the  tools  available  for  the 
study  of  vehicle  crashworthiness.  Although  many  possibilities  sug- 
gest themselves  in  terms  of  further  development  and  refinement 
in  the  analysis  procedure  and  computer  program  (some  of  which 
are  noted  below)  the  most  important  and  potentially  fruitful  ac- 
tivity lies  in  the  immediate  and  continued  application  of  the 
analysis  to  actual  vehicles  and  crash  events.  Such  activity  will 
bring  the  analysis  to  real  vehicle  related  problems  more  quickly 
and  will  serve  as  a guide  to  its  further  development  and  applica- 
tion. Thus  our  primary  recommendation  is  that  provision  be  made 
for  the  continued  and  extensive  use  of  the  program  in  connection 
with  actual  vehicle  behavior.  Some  suggestions  for  such  applica- 
tions as  well  as  for  further  analysis  and  program  development  are 
as  follows: 

Applications 

© Additional  simulation  of  controlled  crash  events 
along  the  lines  of  those  attempted  in  this  work  but 
extended  to  include  machinery  and  drive  train  com- 
ponents and  side- on  impacts . 

© Use  of  the  program  in  connection  with  studies  of 
alternate  vehicle  designs,  including  potential 
structural  modifications  and  configuration  changes. 

© Use  of  the  program  to  generate  structural  crush  data 
for  subsequent  incorporation  in  the  more  simple 
lumped  parameter  or  hybrid  class  of  programs. 

© Comparative  studies  of  the  deformation  and  crash- 
worthiness  of  existing  vehicle  designs. 

Program  Development 

• Development  of  suitable  approximate  moment- curvature 
relations  for  sheet  metal  as  a means  of  improving 
computational  efficiency. 
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• Incorporation  of  plate  and  beam  hinges  at  selected 
locations  in  the  main  version  of  the  program. 

• Development  of  higher  order  and  special  hybrid  element 
formulations  possibly  to  include  tearing  or  wrinkling 
within  the  element. 

• Development  of  an  implicit  rather  than  explicit  in- 
tegration procedure  either  as  an  alternate  to  the 
present  explicit  procedure  or  for  selective  use  for 
part  of  the  structure  or  portions  of  its  response 
history . 
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APPENDIX  A 


TRIANGULAR  PLATE  DISPLACEMENT  AND  CURVATURE  FUNCTIONS 


This  appendix  provides  detailed  expressions  for  the  trans- 
verse displacement  and  curvature  in  the  triangular  plate  elements. 
In  terms  of  triangular  coordinates,  Zienkiewicz  (Ref.  26)  gives 
the  displacement  of  the  plate  as  follows: 


uz  61;1  A B2C2  + 83C3  + 64<-Cl?2  + 2?l52;3-> 

+ ••  • + S9(?l?2  + Z-?li;2?3) 


(Al) 


where  ^ 2.  * ^2  ’ ^3  are  t*ie  triaagular  coordinates  and  8-^  , $2  ’ ‘ ’ ^9 

are  the  constants  depending  on  the  boundary  conditions.  We  define 
the  nodal  rotations  as 
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Since  the  displacements  are  zero  at  the  nodes,  the  unknowns  of 
equation  (Al)  can  be  determined  in  terms  of  nodal  rotations,  and 
equation  (Al)  can  be  expressed  as 
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and  N^, (i=l , . . . , 6)  are  the  shape  functions  given  by  Meek  (Ref.  27). 
Once  the  deformation  is  assumed,  the  slopes  follow 
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and  finally  the  curvatures  are  given  as 
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With  the  element  displacement  defined  by  equation  (Al)  and 
the  slopes  by  equation  (A2),  Meek  (Ref.  27)  gets  the  following 
shape  functions : 

J1  = b2^1C3+2  C1C2C3)  'b3(ClC2  +I  ?1^2C3) 

2 1 2 1 
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When  second  partial  derivatives  of  these  shape  functions  are  ob 
tained,  corresponding  to  the  expressions  in  equation  (A5) , we 
get 
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(A7) (Concl) 


where  A is  the  area  of  the  plate  and  the  constants  are  related  to 
the  coordinates  of  the  node  points  by: 


A A 


A A 


b2 

b3  = 


yryk  ■ 

ai  = Xk"Xj 

A 

> 

a2  = "xk 

A 

-yj  : 

a3  = xj 
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APPENDIX  B 


COMPUTER  PROGRAM  INPUT 

This  appendix  presents  a description  of  the  computer  program 
input  including  definition  of  the  input  quantities,  descriptions 
of  the  user  supplied  driving  subroutine,  FREEFD , and  external 
data  file  options  and  requirements . 

1.  INPUT  FORMAT  FOR  DOT  VEHICLE  STRUCTURES  ANALYSIS  PROGRAM 
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Y(I),Z(I),T(I)  Y(I) , - cross  section  coordinates  of  point  at  begir 

Z(I)  ning  of  segment  I (in.) (See  figure  5) 

T(I)  - segment  thickness  (in.) 
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Repeat  Group  9.1  NPRU  Times;  Omit  for  NPRU 
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2.  EXTERNAL  FILE  REQUIREMENTS 


Restart  unit  - Logical  unit  25 

Permanent  storage  for  checkpoint  data;  approximately  27,000 
words  per  checkpoint  required.  Read  and  write  operations  are  un- 
formatted FORTRAN  type. 

Time  history  unit  - Logical  unit  22 

Permanent  storage  containing  header  records  and  all  output 
quantities  (i.e.,  NPRU,  etc.);  approximately  NPRUxMXSTEP  words  are 
required.  Read  and  write  operations  are  unformatted  FORTRAN  type. 

Mesh  data  file  - Logical  unit  20 

Permanent  storage  containing  initial  mesh  data  (nodes  and 
elements)  and  displacements  for  all  nodes  at  selected  time  steps 
(NPIC) ; the  approximate  number  of  words  required  is 

5 xNUMEL+4x ( NP IC+1 ) xNNODE . 

Read  and  write  operations  are  formatted  FORTRAN  type. 


3.  FORCING  FUNCTION  SPECIFICATION 

Forcing  functions  (i.e.,  displacements  on  nodal  forces  as 
functions  of  time)  are  presently  provided  through  a user  supplied 
subroutine  FREEFD.  The  program  may  be  used  both  for  initializa- 
tion and  subsequent  driving  forces  and  motions  based  on  the  fixed 
point  variable,  KONWAV. 

K0NWAV=2  - initialization 

K0NWAV=1  - forcing  functions 

Displacements,  velocities,  accelerations  and  external  nodal  forces 
are  found  in  UD,  UD1 , UD2  and  FORCD;  counted  six  per  node  in  the 
sequence  u^,  u^ , 9x>  9^_,  9^.  The  program  is  called  at  every 

time  step  and  the  current  valued  time  and  the  time  count  are  con- 
tained in  variables  TIME  and  NTSTEP.  These  may  be  used  as  argu- 
ments in  expressions  for  forces  or  motions  as  functions  of  time. 
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APPENDIX  C 


PROGRAM  DESCRIPTION 
1.  GENERAL 

This  appendix  presents  a description  of  the  computer  program 
developed  in  this  investigation  for  the  dynamic,  large  deflection 
elastic-plastic  response  of  vehicle  structures.  The  program  con- 
tains beam  and  triangular  plate  finite  elements  and  treats  the 
large  deflection  phenomena  by  decomposing  the  element  deforma- 
tions into  rigid  body  and  deformation  components  before  computing 
element  forces.  The  numerical  integration  technique  employed  is 
an  explicit  Newmark  3-method  with  3 = 0 and  a lumped  mass  approxi- 
mation for  the  inertial  properties. 

As  described  in  this  appendix  the  program  consists  of  ap- 
proximately 2000  cards  and  normally  executes  in  65K  words  of  core 
storage  on  the  UNIVAC  1108  computer  under  the  EXEC- 8,  version  27 
operating  system.  A version  of  the  program  has  been  prepared  and 
executed  in  240K  bytes  on  an  IBM  370/158  system  under  the  H 
compiler  with  no  optimization.  The  capacity  of  the  program  as 
described  here  is  approximately  150  nodes,  150  elements,  100  dis- 
placement boundary  conditions,  5 different  sets  of  material  prop- 
erties and  10  distinct  beam  cross  sections.  The  solution  proce- 
dure is  completely  in  core  and  makes  no  use  of  peripheral  devices 
for  temporary  (i.e.,  scratch)  storage  during  execution.  However, 
a restart  capability  is  provided  via  logical  unit  25  and  output 
data  storage  is  provided  via  units  20  and  22.  The  use  of  these 
files  is  optional  and  they  should  be  defined  as  dummy  files  if 
not  required.  The  following  sections  of  this  appendix  provide 
a brief  description  of  the  arrangement  and  function  of  the  various 
subroutines  and  the  allocation  of  the  main  data  storage. 
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2.  SUBROUTINE  DESCRIPTIONS 


The  overall  structure  and  the  names  and  logical  relation 
among  subroutines  is  shown  in  Figure  C-l.  The  subroutines  col- 
lect into  three  main  operating  groups  relating  to  data  input 
QREADIN) , preliminary  calculations  (ASSBLE) , and  the  main  solu- 
tion procedure  (SOLVE)  plus  four  small  subroutines  (utilities) 
which  are  used  throughout  the  program.  Brief  descriptions  of 
the  functions  of  the  individual  programs  follow. 

MAIN.  The  main  program  is  fairly  short  and  contains  di- 
mensioning for  the  main  data  storage  and  calls  to  READIN,  ASSBLE, 
and  SOLVE.  No  computational  or  control  operations  are  carried 
out  in  this  program  except  the  calls  to  these  subroutines. 

READIN . This  program  reads  all  of  the  input  data  from  cards 
(except  for  beam  cross  section  data)  according  to  the  input  for- 
mats described  in  Appendix  B and  initializes  the  data  files  if 
required . 

BGEOM.  Reads  data  describing  the  cross  section  properties 
of  beams  and  generates  certain  additional  data  such  as  segment 
lengths  and  the  torsional  constants. 

ASSBLE . This  program  initilizes  the  element  data  storage, 

assembles  the  mass  matrix  and  the  initial  transformations,  [B  ] 

o 

for  the  nodal  coordinate  systems  and  generates  the  initial  val- 
ues of  the  pointing  vectors  [n]  , [tT]  and  [£]  and  the  nodal  com- 
ponents of  the  rigid  links  [ZT]  . 

TASME . Calculates  contributions  to  the  lumped  mass  matrix 
for  plate  elements  and  forms  the  initial  element  coordinate  sys- 
tem [E]  for  the  plate. 

QUAD . An  algorithm  for  lumping  the  plate  rotational  inertia 
at  the  three  defining  nodes  based  on  the  properties  of  three 
quadrilaterals  formed  in  the  triangle  by  extending  perpendiculars 
from  the  centroid  to  three  sides. 
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Figure  C-l  Computer  Program  Structure 


BASME.  Calculates  contributions  to  the  lumped  mass  matrix 
for  beam  elements  and  forms  the  initial  element  coordinate  system 
[E]  . 

EIGEN.  Generates  the  principal  values  and  directions  for 
the  assembled  rotational  mass  matrix  at  a node.  An  IBM  utility 
program  found  in  Ref.  29. 

SOLVE.  This  program  contains  the  numerical  integration  pro- 
cedure and  the  equations  of  motion  (equation  (2)).  The  principal 
subroutines  called  by  this  program  are  FREEFD  for  external  driving 
forces  or  motions,  FRCIN  for  internal  nodal  forces  and  OUTPUT 
which  provides  printed  output  and  data  files. 

FREEFD.  This  program  provides  nonzero  initial  conditions  or 
external  driving  forces  or  motions  at  specified  nodes  and  is  mod- 
ified for  each  application  as  described  in  Appendix  B. 

RESTRT . This  program  generates  or  reads  checkpoint  data 
files  which  contain  sufficient  information  to  restart  the  solu- 
tion procedure  at  specified  values  of  time  step. 

ROTATE . The  program  is  used  to  constrain  nodal  motions  in 
directions  other  than  parallel  to  the  global  axes . 

FRCIN.  Provides  internal  nodal  forces  for  use  in  the  equa- 
tions of  motion  within  SOLVE.  The  program  updates  the  nodal  co- 
ordinate transformations  [B] , and  calls  subroutines  TFRCIN  for 
plate  forces  and  BFRCIN  for  beam  forces. 

TFRCIN . Calculates  plate  element  deformations  in  the  ele- 
ment coordinate  system,  calculates  nodal  forces  in  the  element 
system  and  transforms  forces  to  the  nodal  coordinate  system. 

VECTOR.  Calculates  the  deformed  length  of  plate  element  side 
and  the  components  of  a vector  normal  to  the  plate  reference  sur- 
face defined  by  the  displaced  positions  of  the  three  node  points. 

CURVAT . Provides  algebraic  expressions  for  curvature  com- 
ponents at  the  midpoints  of  the  three  element  sides  as  functions 
of  the  nodal  rotations . 
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GENFOR.  Controls  the  numerical  integration  of  stresses 
through  the  thickness  of  the  plate  to  obtain  internal  cross  sec- 
tion moments  and  forces. 

MISES . Provides  the  biaxial,  elastic-plastic  stress  strain 
algorithm  described  in  Section  2.4.1. 

BFRC1N.  Calculates  beam  element  deformations  in  the  element 
coordinate  system,  calculates  nodal  forces  in  the  element  system, 
performs  the  operations  associated  with  the  master-slave  rela- 
tions and  transforms  forces  to  the  nodal  coordinate  systems. 

LOCFRC . Calculates  midplane  strains  and  curvatures  for  the 
beam  and  nodal  forces  and  moments  in  the  beam  element  coordinate 
system . 

GENFBM.  Controls  the  numerical  integration  of  stresses  over 
the  cross  section  of  the  beam  to  obtain  internal  cross  section 
moments  and  forces. 

STRES . Provides  the  uniaxial,  elastic-plastic  stress  strain 
algorithm. 

OUTPUT . Provides  printed  output  of  those  quantities  selected 
for  output,  and  updates  the  time-history  and  deformed  geometry 
data  files. 
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3.  DATA  ALLOCATION 


The  allocation  of  the  main  items  of  data  throughout  the 
program  is  as  follows. 


COMMON  (blank) 

XC (300) 

YC (300) 

ZC (300) 

UD (900) 

UD1 (900) 

UD2 (900) 
FORCD (900) 

E (16, 5) 

IX (9 , 300) 

NODDIS (100) 

NUMNP 

NUMEL 

NDGREE 

MUD 

NUMDIS 

COMMON /CONTRL/ 
KONTRL(IO) 

COMMON/ DYNAM/ 
DELT 
TIMEND 
MXSTEP 
NTSTEP 
TIME 


--  Undeformed  nodal  coordinates 


--  Nodal  displacement,  velocity  and  acceler- 
ation (six  per  node) 

--  External  nodal  forces 

--  Properties  for  up  to  five  different 
materials 

--  Array  of  element  node  correspondences  and 
other  element  parameters 

--  Array  of  packed  word  constraint  data 

--  Number  of  node  points 

--  Number  of  elements 

--  Degrees  of  freedom  per  node  (6) 

--  Not  used 

--  Number  of  displacement  points 

--  Control  flags 

--  Time  increment 
--  Maximum  time  value 
--  Maximum  number  of  time  steps 
--  Current  step  count 
--  Current  value  of  time 


The  names  of  arrays,  and  number  of  subscripts  differ  in  some 
subroutines  with  correspondence  arranged  through  subroutine 
calling  arguments  or  explicit  EQUIVALENCE  statements. 
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COMMON/ XSECT/ 
NSEG(IO) 

Y (20 , 10) 

Z (20 , 10) 

. XLEN(20, 10) 
T(20 , 10) 
TORCON(IO) 
KLOS (10) 
ABM(IO) 
FIXX(IO) 
FIYY(IO) 
FIZZ (10) 


Number  of  segments  in  beam  cross  sections 
Beam  cross  section  coordinates 

Beam  cross  section  segment  lengths 
and  thicknesses 

Torsional  constant 
Closed  wall  section  flag 

Cross  section  area  and  moments  of  inertia 
about  x,  y and  z axes 


COMMON  /STR/ 

II 

INDEX (300) 

STRAIN(45) 
STRESS (45) 
SIDE(900) 

IP  (6) 

STRS (10000) 


--  Number  of  plate  cross  section  integration 
points 

--  Array  of  pointers  to  first  item  in  element 
stress  array  STRS 

--  Working  storage  for  stress  and  strain 
calculations 


--  Initial  lengths  of  plate  element  sides 
(three  per  plate  element) 

--  Plate  element  side/node  counter 


--  Element  stress  and  strain  storage-  Format 
for  plates  is  3x11x7  with  each  block  of  7 
laid  out  as 


Te  , £ , £ , o , o , o , o 

1 xx ’ yy  xy  xx ’ yy  xy  y 


for  II  points  on  3 sides  for  each  element. 
Format  for  beams  is  2xNPTSx4  with  each 
block  of  4 laid  out  as 


lexx’°xx' Vblank] 

for  NPTS  points  at  2 ends  of  the  beam.  Note 
that  NPTS=NSEG+KLOS. 


COMMON/ JUNK/ 

TL(9 , 300)  Components  of  the  [B]  transformation  laid 

out  by  columns  of  three,  three  columns 
per  node 

POINT (3,3,300)- * _C ompo n en t s of  initial  plate  normal  in 

[x,y,z]  system 
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COMMON/ BEAMS/ 

RMASS (3 ,3,100)  - 

DICOS(3, 3, 150)  - 

BAR (12 00) 


PMASS (4 , 100) 
NSECT 

COMMON/ OUTPA/ 

NPRU 

NPRS 

NPFREQ 

NPIC 

OUT (9,40) 
NPOUT (2,10) 
TITLE (20) 


- Components  of  assembled  rotational 
inertia  tensor  before  diagonal ization 

- Components  of  initial  element  trans- 
formation [E] 

- Beam  and  rigid  link  storage  laid  out  as 
[^o’^o;AG;^  = 12  words 

at  two  ends  of  beam  = 24  words  per  beam 

- Added  nonstructural  lumped  mass  at  nodes 
one  translational  and  three  rotational 
components  per  node 

- Number  of  beam  cross  section  types 


- Number  of  displacement  output  points 

- Number  of  stress  output  points 

- Print  frequency 

- Number  of  picture  output  intervals 

- Packed  word  output  list 

- Picture  output  list 

- Problem  title 
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UCK2)s?(SCn*C(3)+CU)*SC2)*S(3))*UC(mCU)*C(2)*UCe2)  00003300 
U<S(1)*$(3)»C<U*S(2)*C(J))*UC(J)  0 0003400 
U(K3)=.C(2)*SC3)*UCC1)*S(2)*UC(2)*CC2)*C(3)*UC(3)  00003500 
RETURN  00003600 
END  00003700 
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116.  DA»  ABS  ( A 1 1 + A 1 ( I ) ) OOOU0OO 

119.  IF  ( DA  .IE,  0.0)  GO  TO  425  00011900 

120.  SUM3SUM+ABS(A1I*A1 (I))/DA  00012000 

121.  AlCDsAlI  00012100 
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1.  INTRODUCTION 


The  increased  concern  in  recent  years  over  vehicle  safety  has 
focused,  in  part,  on  the  ability  of  a vehicle  to  sustain  a crash 
event;  to  absorb,  redirect  and  otherwise  manage  the  severe  forces 
and  energy  associated  with  the  event;  and  thereby,  to  lessen  the 
severity  of  the  environment  to  which  passengers  are  exposed.  In 
addition  to  basic  vehicle  crashworthiness,  serious  attention  has 
also  been  devoted  to  such  factors  as  the  propensity  of  a vehicle 
to  damage  in  minor  crash  events  and  the  threat  posed  to  pedestrians 
or  other  vehicles  by  particular  vehicle  designs. 

Since  quantitative  knowledge  of  the  deformation  characteris- 
tics of  vehicle  structures  is  an  intrinsic  requirement  in  all  con- 
siderations of  this  type,  it  is  not  surprising  that  the  developing 
concern  for  vehicle  safety  has  given  rise  to  intense  and  continuing 
efforts  to  provide  mathematical  descriptions  of  structural  deforma- 
tions in  the  crash  environment.  This  report  presents  the  results 
of  a research  project  undertaken  by  IIT  Research  Institute  (IITRI) 
for  the  Department  of  Transportation  to  develop  a finite  element 
computer  program  for  use  in  the  dynamic  analysis  of  vehicle  struc- 
ture, including  sheet  metal,  in  a crash  environment. 

Present  analytical  procedures  fall  largely  into  two  broad 
categories  which  may  be  termed  lumped  parameter  or  equivalent  sys- 
tems and  formal  approximation  techniques  (primarily  finite  element). 
In  the  former  category,  the  vehicle  is  treated  as  an  assemblage 
of  lumped  masses  and  resistance  elements  which  are  selected  to 
represent  specific  subassemblies  and  components.  The  properties 
of  the  various  elements  of  such  models  are  determined  by  a static 
and  dynamic  crush  testing,  direct  measurement,  and  the  judicious 
use  of  simple  analysis  procedures.  Such  techniques  are  valuable 
in  that  they  provide  an  efficient  and  inexpensive  method  of  simu- 
lation, are  fairly  accurate  for  well  established  vehicle  configura- 
tions and,  in  effect,  serve  as  a depository  for  previous  test  and 
design  results.  The  limitations  of  this  type  of  procedure  are 
principally  that  they  are  heavily  reliant  on  previous  testing; 
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large ly  limited  to  well-established  vehicle  configurations  and 
only  weakly  tied  to  basic  principles  of  structural  behavior. 

The  finite  element  approach  represents  an  attempt  to  make 
up  for  the  limitations  inherent  in  the  lumped  parameter  analyses, 
by  employing  more  formal  approximation  techniques  in  the  discre- 
tization of  the  structure  and  a greater  reliance  on  more  funda- 
mental structural  representations  and  properties.  The  limitations 
of  this  approach  are  found  in  the  inherent  tendency  toward  more 
complicated  and  expensive  computations  and  the  difficulties  found 
in  modeling  the  extensively  complex  phenomena  associated  with  the 
crash  environment.  Such  phenomena  include  large  deflections  and 
rotations  in  the  deformed  structure,  regions  of  intense  curvature 
(wrinkling),  material  strain  rate  effects  and  interference  and  con- 
tact among  structural  components  during  the  response.  These  pro- 
cedures, indeed,  are  not  totally  free  of  a reliance  on  testing  and 
analytical  judgment  and,  in  fact,  can  be  viewed  as  a somewhat  more 
formal  lumped  parameter  system  and  used  in  connection  with  the  more 
empirical  procedure  in  hybrid,  combination  models. 
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2.  ANALYSIS 


The  response  of  vehicle  structures  under  crash  loadings  is  a 
complex  process  primarily  involving: 

• transient,  dynamic  behavior 

• complicated  framework  and  shell  assemblages 

• large  deflections  and  rotations 

• extensive  plastic  deformation 

Previous  attempts  at  a formal  analysis  of  this  process  have  been 
only  partly  successful  due  to  a variety  of  limitations  which,  in 
particular  instances , have  included  inadequacies  in  element  formu- 
lations , material  representations  or  solution  procedures.  The 
work  presented  in  this  report  represents  an  attempt  to  develop  a 
finite  element  program  which  is  specifically  tailored  to  the  class 
of  problems  inherent  in  vehicle  crash  response,  and  which  employs 
or  extends  current  avenues  in  finite  element  analysis  which  seem 
best  suited  to  such  problems. 

In  this  section  major  features  of  the  present  work  are  brief- 
ly described,  and  some  rationale  is  offered  for  their  use  in  the 
context  of  vehicle  analysis. 

2 . 1 Solution  Procedure 

The  analysis  procedure  is  based  on  explicit,  timewise  numer- 
ical integration  of  the  equations  of  motion  for  the  node  points 
(three  translations  and  three  rotations  per  node).  The  choice  of 
an  explicit,  rather  than  implicit,  procedure  is  exploited  through- 
out the  analysis  by  using  "lumped"  nodal  masses  and  by  calculating 
the  internal  forces  at  the  nodes  from  direct  integration  of  the 
element  stress  fields  without  reference  to  element  or  assembled 
stiffness  matrices  for  the  structure. 

The  equations  of  motion  at  a typical  node  point  for  transla- 
tions and  rota tions,  respective ly , are  written  as  follows: 

Trans lations 

[M][ug]  = [fg]  - [fG] 
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(1) 


where 


[M]  - the  diagonal,  lumped  mass  matrix  for  the  node 

[ug]  fuxG’ uyG  ,uzG^  is  tbe  acceberat:i-on  vector  for 

the  node  referred  to  a global  system  [Xn , Y^, , Zn  ] 
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[f^]  - the  vector  of  external,  applied  loads  at  the 
node  referred  to  the  global  system 
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wnere  all  components  are  referred  to  a coordinate  system  [x,y,z] 
wi'iich  coincides  with  the  principal  axes  of  the  lumped  mass  moments 
of  inertia  at  each  node  and  which  rotates  with  the  node  (i.e.,  a 
set  of  rigid  body  axes  for  each  node)  and  where  [I  ,1,1]  are  the 

X y z 

principal  mass  moments,  [a  , a , a ] are  instantaneous  angular  ac- 
x y z 

celerations,  [oj  ,oo  , w ] are  the  angular  velocities,  and  where  [m  ] 

x y z 

and  [m  ] are  the  external  and  internal  moment  vectors  at  the  node 
point . 


The  orientation  of  the  nodal  axes  [x,y,z]  with  respect  to  the 
global  axes  at  any  time  during  the  motion  is  established  by  the 
components  of  three  unit  vectors  b-^,b2,b^  which  remain  fixed  along 
the  nodal  axes  [x,y,z],  respectively,  as  the  node  rotates.  If  the 
global  components  of  these  three  unit  vectors  are 

C:  (b  ~>b  --b  ~) 

1 xx’  yx  zx' 

b 0 : (b  — , b — , b — ) 

2 xy  yy  zy 

bQ:  (b  — ,b  — ,b  — ) 
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(3a) 

(3b) 

(3c) 


-i. 

then  the  components  of  any  vector,  v,  transform  from  the  nodal  to 
global  system  by  the  transformation 

[vG]  = [B][v]  (4) 

where  the  columns  of  [B]  are  the  components  of  the  three  unit  vectors. 

The  choice  of  an  explicit  integration  procedure  and  a direct 
calculation  of  nodal  forces  results  in  a program  with  minimum  (but 
still  substantial)  computer  storage  requirements.  This,  in  turn, 
is  equivalent  to  a capability  of  processing  reasonably  detailed  and 
extensive  structural  models  with  relative  ease.  An  early  capability 
for  handling  relatively  large  problems  is  considered  essential  in 
developing  realistic  simulations  of  actual  crash  events  because  of 
the  complex  geometry  and  construction  found  in  real  vehicles.  Fur- 
thermore, this  procedure  lends  itself  to  a reasonably  well  organized 
modular  program  which  admits  subsequent  change  and  development,  in- 
cluding the  further  incorporation  of  an  implicit  method. 

2.2  Element  Formulation 

The  treatment  of  large  displacements  and  rotations  employs  a 
decomposition  of  the  element  displacement  field  into  a rigid  body 
rotation  and  translation  associated  with  a local  coordinate  system 
attached  to  and  moving  with  the  element,  and  a remaining  displace- 
ment field  which  describes  the  deformation  of  the  element  relative 
to  the  current  position  of  the  element  axes.  This  transformation, 
in  effect,  removes  the  average  rigid  body  rotation  of  the  element 
and  allows  the  use  of  small  or  moderate  deflection  element  formu- 
lations in  the  calculation  of  element  forces.  In  this  manner, 
extremely  large  rotations  and  deflections  can  be  accommodated  by 
the  analysis  with  accuracy  depending  primarily  on  the  size  of  the 
elements  relative  to  the  curvature  of  the  structure.  Although  a 
formal  convergence  theorem  is  lacking,  the  decomposition  does  hold 
for  infinitesimal  regions  and  numerical  studies  show  excellent 
agreement  with  classical  solutions.  The  computer  program  at  pres- 
ent includes  low  order  triangular  plate  elements  and  three- 
dimensional  beam  elements. 
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A coordinate  system  [x,y,z]  is  embedded  in  each  element  and 
serves  to  define  the  rigid  body  rotation  of  the  element  and  as  a 
reference  for  element  distortions  and  forces.  The  components  of 
a vector,  V,  transform  from  the  element  coordinate  system  to  the 
global  coordinate  system  by  the  time  varying  transformation 

[vB]  = [ E ] [ v] 

where  the  elements  of  E are 
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xy 
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which  correspond  columnwise  with  the  global  components  of  unit  vec- 
tors e-^,e2,e2  oriented  with  the  element  axes  [x,y,z]  respectively. 
An  element  coordinate  system  for  a triangular  plate  in  the  unde- 
formed and  deformed  positions  is  shown  in  Figure  1.  The  transfor- 
mations [E]  for  each  element  and  [B]  for  each  node  together  with 
the  components  of  the  normal  vector,  n,  at  each  node  or  an  ele- 
ment, provide  sufficient  information  to  calculate  deformations 
in  the  element.  The  deformations,  in  turn,  are  used  to  calculate 
stresses  within  the  element  and,  subsequently,  forces  and  moments 
at  the  nodes.  Finally,  the  nodal  forces  and  moments  are  entered 
in  the  appropriate  summations  in  the  equations  of  motion  (equa- 
tions ( 1)  and  (2)  ) . 


2.3  Material  Properties 

The  computer  program  currently  uses  simple  elastic-plastic 
stress  strain  laws;  a uniaxial  relation  for  beam  elements  and  a 
biaxial  strain  hardening  Mises  model  for  plates.  Element  forces 
and  bending  moments  for  given  strain  fields  are  calculated  by 
piecewise  linear  numerical  integration  of  the  stresses  at  selected 
points  in  the  cross  section.  The  program  is  designed  to  accommo- 
date other  constitutive  relations  readily  or  to  accept  explicit 
relations  among  thrust,  moment,  extension  and  curvature.  In  fact 
linear  relationships  are  provided  at  this  level  which  result  in 
more  efficient  computations  for  this  class  of  problem. 
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Figure  1 Three-Dimensional  Plate  Element 
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3.  ANALYSIS  RESULTS 


A substantial  number  of  test  and  demonstration  problems  were 
treated  during  the  development  of  the  analysis  and  computer  pro- 
gram reported  here.  For  the  most  part  the  test  problems  dealt 
with  relatively  simple  structural  configurations  having  known  re- 
sults and  were  selected  to  isolate  and  test  particular  aspects  of 
the  formulation  and  program  during  development.  Thirteen  such 
test  problems  of  successively  increasing  complexity  were  formally 
reported  in  progress  reports . In  all  cases  the  analyses  were  pur- 
sued until  satisfactory  comparison  with  previous  results  were  ob- 
tained or  at  least  to  the  point  at  which  differences  between  re- 
sults could  clearly  be  reconciled  on  the  grounds  of  differences  in 
loadings,  boundary  conditions,  or  ambiguities  in  the  results  used 
for  comparison  purposes . The  contents  of  these  test  problems  are 
recapitulated  briefly  in  Table  1. 

In  addition  to  the  more  idealized  test  problems  an  effort  was 
made  using  relatively  large  scale  simulations  of  actual  vehicle 
crash  events.  The  events  which  are  the  objects  of  these  simula- 
tions were  selected  from  a series  of  barrier  tests  conducted  by 
Dynamic  Science  which  involved  frontal  impacts  at  various  velocities 
of  a 1968  Plymouth  Fury  having  various  portions  of  the  front  end 
structure  and  machinery  removed  in  different  tests.  The  object  of 
these  tests  was  to  isolate  and  identify  the  dynamic  resistance 
properties  of  the  various  front  end  components  in  three  groups, 
namely,  frame  structure,  sheet  metal  structure  and  the  engine  and 
machinery  components.  Simulations  were  developed  for  two  test 
configurations  involving  the  stub  frame  alone  (test  283-51)  and  the 
stub  frame  and  sheet  metal  together  (test  283-53) . 

The  finite  element  models  developed  for  these  two  test  con- 
figurations represent  one-half  of  the  structural  components  for- 
ward of  the  vehicle  fire  wall  and  are  based  on  an  assumption  of 
symmetry  for  both  the  vehicle  and  its  response  modes  about  the 
center  plane. 
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SUMMARY  OF  TEST  PROBLEMS 
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The  frame  model  consists  entirely  of  beam  elements  and  the  frame 
and  sheet  metal  model  contains  the  frame  model,  plate  elements 
representing  the  sheet  metal,  and  additional  beam  elements  repre- 
senting stiffener  sections  in  the  hood.  In  both  cases  the  appro- 
priate inertial  properties  of  the  remainder  of  the  vehicle  are 
provided  as  added  masses  at  the  vehicle  center  of  gravity  and  are 
attached  to  the  frontal  structure  by  means  of  rigid  links. 

As  an  illustration,  we  briefly  review  the  results  obtained  from 
a simulation  of  the  stub  frame  and  sheet  metal  barrier  test  (test 
233-53;  44  mph) . An  isometric  view  of  the  finite  element  mesh 
for  the  stub  frame  and  sheet  metal  model  before  impact  is  shown  in 
Figure  2.  This  model  comprises  the  left  front  quarter  of  the 
1968  Plymouth  employed  in  the  test.  The  major  components  of  the 
model  are  the  hood,  left  front  fender,  wheel  well  and  frame  and 
radiator  support.  The  wheel  well  is  shown  separately  for  the  sake 
of  clarity  in  the  figure  but  is  actually  an  integral  part  of  the 
model  and  is  attached  along  the  lines  defined  by  the  match  points 
1 through  4 indicated  in  the  figure.  The  thickness  of  all  sheet 
metal  elements  was  taken  as  0.039  in.  material.  The  hood  portion 
of  the  model  is  free  to  move  relative  to  the  fender  at  their  line 
of  intersection,  and,  in  addition  to  the  plate  elements  shown, 
contains  2.0x1.0x0.05  in.  hat  section  stiffeners  at  reinforce- 
ments. As  with  the  stub  frame  model  additional  mass  was  lumped  at 
the  approximate  location  of  the  wheel  suspension  and  also  at  the 
vehicle  center  of  gravity.  In  this  case  all  points  along  the  rear 
edge  of  the  fender  and  hood  as  well  as  the  end  of  the  frame  were 
attached  to  the  vehicle  center  of  gravity  by  means  of  rigid  links. 
The  model  was  driven  by  an  initial  longitudinal  velocity  of  774  in./ 
sec  for  all  nodes  with  longitudinal  constraints  applied  at  nodes  in 
the  forward  section  of  the  model  successively  as  they  came  into 
contact  with  the  barrier. 

The  deformed  shape  of  the  model  at  20  msec  after  impact  is 
shown  in  Figure  3.  At  this  point  the  forward  motion  of  the  center 
of  gravity  is  approximately  15  in.  and  the  forward  portions  of  the 
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frame  assembly  and  the  first  two  lines  of  nodes  in  the  hood  are  in- 
dicated in  the  forward  portions  of  the  model  with  the  bumper  and 
bumper  brackets  collapsed  downward  and  the  hood  and  fender  collapsing 
upward  and  outward  respectively.  Separation  at  the  hood  and  fender 
junction  line  is  clearly  evident  and  the  rear  portion  of  the  hood 
has  begun  to  buckle  upward  away  from  the  fender.  The  wheel  well 
is  not  greatly  involved  although  the  forward  wall  has  bent  to  con- 
form to  the  distortions  in  the  forward  section  of  the  fender. 

Longitudinal  acceleration  at  the  center  of  gravity  and  total 
barrier  force  as  obtained  from  test  and  analysis  are  shown  in  Fig- 
ures 4 and  5.  In  both  cases  the  test  displays  the  effects  of  filter- 
ing which  somewhat  obscures  the  comparisons.  The  calculated  accel- 
erations show  an  initial  peak  of  about  10  g associated  with  collapse 
of  the  bumper  and  its  supports  followed  by  a larger  peak  of  30  g 
as  the  fender  and  hood  are  involved.  Maximum  acceleration  in  the 
test  data  is  about  25  g.  The  calculated  barrier  force  shows  an  ini- 
tial plateau  at  about  40,000  lbs  followed  by  a sudden  rise  to  about 
80,000  lbs  as  the  hood  and  fender  make  contact  with  the  barrier. 

The  corresponding  maximum  barrier  load  from  test  (50  percent  of 
full  value)  is  about  60,000  lbs  during  this  phase  of  the  motion. 

The  test  simulations  presented  here  are  by  no  means  a complete 
or  exhaustive  study  even  for  the  limited  series  of  tests  which  were 
used  for  comparison  purposes.  For  a variety  of  reasons,  the  re- 
sults given  apply  only  to  the  initial  phases  of  the  respective 
crash  events  and  improvements  could  also  be  made  in  the  models  (the 
provision  of  rear  wheels  or  other  vertical  constraints  at  the  cen- 
ter of  gravity,  for  instance)  and  in  the  associated  data  processing. 
Despite  these  limitations,  however,  the  results  appear  to  be  a rep- 
resentative portrayal  of  the  phenomena  associated  with  vehicle  bar- 
rier collisions  generally,  and  a reasonable  correlation  with  the 
particular  test  data  for  which  comparisons  were  attempted. 
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Figure  2 Stub  Frame  and  Sheet  Metal  Model  (Undeformed) 


© 


15 


Figure  3 Stub  Frame  and  Sheet  Metal  Model  (Deformed,  20  msec) 
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Figure  4 Stub  Frame  and  Sheet  Metal  Vehicle  Center  of  Gravity  Acceleration 
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Figure  5 Stub  Frame  and  Sheet  Metal  Total  Barrier  Force 


4.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  report  presents  the  results  of  a research  project  under- 
taken by  IITRI  for  DOT  to  develop  a finite  element  computer  program 
for  use  in  the  dynamic  analysis  of  vehicle  structures,  including 
sheet  metal,  in  a crash  environment.  This  research  program  in- 
volved the  following  major  tasks: 

• A technique  was  developed  for  the  finite  element  analy- 
sis of  the  dynamic  response  of  plate  beam  structures 
involving  very  large  displacements  and  rotations  and 
elastic-plastic  material  behavior.  The  principal  fea- 
ture of  this  technique  involves  the  decomposition  of 
the  element,  displacement  field  into  rigid  body  compo- 
nents and  deformation  components  thus  allowing  the  use 

of  a small  deflection  element  formulation  in  the  analysis. 

« A computer  program  was  developed  incorporating  this  an- 
alysis procedure  for  beam  and  plate  elements  and  rigid 
links  together  with  appropriate  time  integration  pro- 
cedures and  material  property  descriptions. 

# A substantial  number  of  test  and  demonstration  problems 
were  analyzed  with  this  computer  program  ranging  from 
simple  classical  solutions  for  beams  and  plates  through 
large  scale  simulations  of  actual  crash  tests. 


4. 1 Conclusions 

Based  on  the  work  carried  out  in  this  research  program  our 
principal  conclusions  are  as  follows: 

(1)  The  analysis  formulation,  while  somewhat  novel,  is  a for- 
mally sound  procedure  and  a proper  and  useful  approximation  tech- 
nique for  the  dynamic  analysis  of  beam  and  plate  structures 
involving  large  deflections  and  rotations. 

(2)  The  computer  program  developed  is  a correct  rendering  of 
this  analysis  technique  and  provided  accurate  results  with  reason- 
able efficiency  in  comparison  to  available  known  solutions. 

(3)  Finally,  and  most  importantly,  the  computer  program  is 
readily  applicable  to  realistic  vehicle  structures  and  provides 
credible  results  for  actual  crash  events  as  demonstrated  by  the 
simulations  of  the  end-on  barrier  tests. 
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4.2  Re  commendat ions 


In  our  opinion  the  analysis  and  computer  program  resulting 
from  this  research  is  a substantial  improvement  over  previous  at- 
tempts at  a finite  element  analysis  of  vehicle  structures  and  is 
a potentially  valuable  addition  to  the  tools  available  for  the 
study  of  vehicle  crashworthiness.  Although  many  possibilities 
suggest  themselves  in  terms  of  further  development  and  refinement 
in  the  analysis  procedure  and  computer  program  (some  of  which  are 
noted  below)  the  most  important  and  potentially  fruitful  activity 
lies  in  the  immediate  and  continued  application  of  the  analysis 
to  actual  vehicles  and  crash  events.  Such  activity  will  bring 
the  analysis  to  real  vehicle  related  problems  more  quickly  and 
will  serve  as  a guide  to  its  further  development  and  application. 
Thus  our  primary  recommendation  is  that  provision  be  made  for  the 
continued  and  extensive  use  of  the  program  in  connection  with  ac- 
tual vehicle  behavior.  Some  suggestions  for  such  applications  as 
well  as  for  further  analysis  and  program  development  are:: 

App lications 

« Additional  simulation  of  controlled  crash  events 
along  the  lines  of  those  attempted  in  this  work  but 
extended  to  include  machinery  and  drive  train  com- 
ponents and  side -on  impacts. 

© Use  of  the  program  in  connection  with  studies  of 
alternate  vehicle  designs,  including  potential 
structural  modifications  and  configuration  changes. 

• Use  of  the  program  to  generate  structural  crush  data 
for  subsequent  incorporation  in  the  more  simple 
lumped  parameter  or  hybrid  class  of  programs. 

« Comparative  studies  of  the  deformation  and  crash- 
worthiness  of  existing  vehicle  designs. 

Program  Development 

® Development  of  suitable  approximate  moment- curvature 
relations  for  sheet  metal  as  a means  of  improving 
computational  efficiency. 

© Incorporation  of  plate  and  beam  hinges  at  selected 
locations  in  the  main  version  of  the  program. 
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• Development  of  higher  order  and  special  hybrid  element 
formulations  possibly  to  include  tearing  or  wrinkling 
within  the  element. 

• Development  of  an  implicit  rather  than  explicit  in- 
tegration procedure  either  as  an  alternate  to  the 
present  explicit  procedure  or  for  selective  use  for 
part  of  the  structure  or  portions  of  its  response 
history . 
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